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A  FINITE-DIFFERENCE  PROGRAM  FOR  STRESSES  IN 
ANISOTROPIC,  LAYERED  PLATES  IN  BENDING 

INTRODUCTION 


One  critical  feature  associated  with  structural  composites  of  laminated 
construction,  using  materials  or  geometrical  arrangements  that  exhibit  different  elastic 
properties  from  layer  to  layer,  is  the  possibihty  that  the  glued  layers  will  separate  or 
delaminate.  This  was  undoubtedly  realized  from  the  outset  of  their  use,  and  a  brief 
historical  sketch  of  the  American  scene  is  presented  by  Pipes  [  1  ] .  However,  the  earliest 
serious  investigation  into  the  cause  of  delamination-type  failure,  namely  the  interlaminar 
stress  problem,  was  apparently  done  in  Japan  by  Hayashi  [2,3],  who  reported  that  the 
maximum  interlaminar  shearing  stresses  occurred  at  the  free  edge  of  a  laminate  under 
tension.  Hayashi  used  a  plane  stress  model  for  the  layers  and  approximated  the 
interlaminar  shears  by  a  strain-averaging  technique.  Using  a  similar  model,  Puppo  and 
Evensen  [4]  likewise  discovered  a  sharp  rise  in  the  interlaminar  stresses  near  a  free  edge. 
Notably,  the  use  of  the  above  models  ignored  the  interlaminar  normal  stress.  In  two 
publications.  Pipes  and  Pagano  [5,6]  developed  a  finite-difference  program  to  solve  the 
exact  elasticity  equations  for  a  long  laminate  in  uniaxial  extension.  In  their  development, 
the  stresses  are  assumed  independent  of  the  axial  coordinate  and  include  all  six 
components.  The  results  of  this  investigation  show  that  a  sharp  rise  in  both  the 
interlaminar  shear  stresses  and  the  normal  stress  occurs  near  the  free  edge.  Thereafter, 
Oplinger  [7]  did  an  analysis  of  angle  ply  laminates  in  tension  using  a  model  similar  to 
that  of  References  2  through  4.  His  approach  allows  a  large  number  of  layers  to  be 
considered.  Indeed  it  was  discovered  that  a  singularity  in  the  interlaminar  shear  occurs  at 
the  free  edge  of  a  laminate  of  one  particular  type  of  construction.  An  alternative  solution 
to  that  employed  in  the  above  references  is  used  by  Rybicki  [8]  who  applied  a 
three-dimensional  finite  element  formulation.  His  results  agree  with  References  5  and  6. 

The  present  report  marks  the  initial  phase  of  a  study  of  the  interlaminar  stresses 
induced  in  a  layered  laminate  by  bending.  Following  the  approach  used  by  Pipes  [5],  the 
laminate  is  modeled  as  a  continuum  and  the  resulting  elasticity  equations  are  solved  using 
the  finite-difference  method.  This  solution  technique  is  made  possible  by  assuming  that 
the  laminate  is  bent  into  a  cylindrical  surface  such  that  the  stresses  are  independent  of 
the  axial  coordinate.  The  objective  of  this  report  is  to  set  forth  the  mathematical 
framework,  present  some  preliminary  results,  and  to  avail  the  computer  program  to 
others.  The  results  reveal  a  simphfying  symmetry  relationship  in  the  displacements  that 
will  allow  significant  reduction  in  the  size  of  certain  numerical  problems.  In  addition, 
trends  in  the  interlaminar  stress  distribution  are  somewhat  similar  to  those  found  for 
stretching  problems,  in  that  a  sharp  rise  occurs  at  the  free  edge. 


PROBLEM  FORMULATION 


Laminate  Description 

The  laminated  composites  considered  in  this  report  consist  of  rectangular  laminae 
symmetrically  stacked  with  respect  to  a  midplane  and  bonded  together  to  form  a  flat 
laminate.  The  bonding  is  assumed  to  provide  perfect  adhesion  between  the  laminae, 
which  nullifies  the  possibility  of  slip  between  adjacent  laminae  thus  establishing  the 
conditions  of  continuous  displacements  and  tractions  at  each  interface.  Each  individual 
lamina  is  considered  to  be  elastic,  homogeneous,  and  orthotropic  (i.e.,  each  lamina 
possesses  three  planes  of  elastic  symmetry).  The  assumption  of  homogeneity  eliminates 
micro  mechanical  effects  such  as  those  involving  fibers  or  matrix.  The  geometry  of  a 
typical  lamina  and  laminate  is  illustrated  in  Figure  1 .  One  may  note  that  the  orthotropic 
coordinate  axes  (x',y',z)  of  a  lamina  are  referred  through  a  clockwise  rotation  about  z  to 
the  fixed  coordinate  axes  (x,  y,  z)  of  the  laminate.  The  laminae  are  stacked  along  z  to 
form  a  laminate  whose  sides  are  normal  to  x,  y,  and  z.  Each  lamina  is  given  a  layer 
number  m. 

Limiting  the  analysis  to  linear  elastic  materials,  the  constitutive  relation  for  each 
lamina  referred  to  the  x,  y,  z  coordinate  system  is 


where  the  elastic  constants  cy  are  related  to  the  nine  orthotropic  constants  cy  through 

the  well  known  transformation  equations  of  References  9  and  10.*  By  associating  the 
displacements  u,  v,  and  w  with  x,  y,  and  z,  respectively,  the  strains  for  each  lamina  are 
defined  as 

1.  In  using  the  transformation  equations  in  References  9  and  10  substitute  -6  for+0  since 
here  the  constants  are  referred  to  the  unprimed  coordinate  axes  of  the  laminate. 
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where  the  comma  denotes  partial  differentiation. 


Loading  and  Field  Quantities 

Consider  a  laminate  loaded  by  bending  about  y  at  the  ends  x  =  constant. 
Assuming  that  the  laminate  is  long  enough  in  the  x-direction  and  that  Saint-Venanfs 
principle  holds  for  a  laminate,  the  resulting  stress  distribution  will  be  independent  of  x  in 
regions  sufficiently  removed  from  the  areas  of  loading.  Using  this  assumption  and 
following  Lekhnitskii  [11],  the  elastic  strain-stress  relations  can  be  integrated  to  yield 
displacements  for  each  lamina  of  the  form 


U*”  =  (CiY  +  CjZ  +  C3)x  +  U”^(y,  z) 

V*”  =  -■!  +  C4XZ  +  V^(y,  z) 

w™  =  -|  Cjx^  -  C4xy  +  W'T’Cy.z)  ,  (3) 

where  U^,  V^,  and  are  unknown  functions  of  y,  z.  The  layer  number,  m,  is  left  off 
the  constants  Cj  because  it  results  that  each  Cj  must  be  the  same  for  every  lamina  in 

order  to  satisfy  the  displacement  continuity  conditions  at  the  interfaces.  Thus,  the  Cj  are 

found  to  be  properties  of  the  entire  laminate.  The  displacement  equations  (3)  represent 
the  full  three-dimensional  elasticity  solution  that  holds  for  all  points  in  the  laminate. 

To  evaluate  the  Cj,  the  scheme  is  as  follows.  Since  equations  (3)  hold  for  all 

points  in  the  laminate,  they  must  converge  to  the  plane  stress  solution,  which  is  an  exact 
solution,  in  the  interior  region  of  the  laminate.  Integrating  the  relation  [10,12] 


Ci  =  B[jMj 


zDjjMj 


ij  =  1>  2,  6 


(4a) 


3 


for  the  case  where  Mi  =  -M  and  Mj  =  Mg  =  0,  the  plane  stress  displacements  are  found  to 
be 

=  (-D'liMz  -  BiiM)x  -  B^iMy  -  ^D'leMyz  +  f(z) 

2 

Vpg  =  -^DifiMxz  -  (B21M  +  Di2Mz)y  +  g(z) 

Wps  =  ^D'liMx^  +  ^D'lgMxy  +  |^Di2My^  +  f*(x)  +  g*(y)  ,  (4b) 

where  By  and  Dy  are  laminate  properties  defined  in  Appendix  A,  and  M  is  the  applied 
moment.  Comparing  equations  (3)  and  (4b)  leads  to  the  results: 

C,  =  0  Cj  =  -D'nM 

(5) 

C3  =  -B'liM  C4  = 

and 

U”^(y,  z)  ^  B^y  +  C4yz  +  U*^(y,  z) 

V™(y,  z)  ->  B^y  +  D^yz  +  V"i(y,  z) 

W™(y,  z)  ->■  -  ]-  D^y^  +  W^(y,  z)  ,  (6) 

where  ^ 

B^  =  -B;,M  ,  By  =  -B^iM  ,  and  Dy  =  -D;2M  .  (7) 

2.  The  extended  forms  (6)  for  U^,  V™,  and  are  not  necessary  to  the  solution. 
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Substituting  the  results  (6)  into  equations  (3)  yields  displacements  of  the  following 
functional  form  for  each  layer 

u”’  =  (CjZ  +  C3)x  +  (By  +  C4z)y  +  U’^Cy,  z) 

=  C4XZ  +  (By  +  DyZ)y  +  V”^(y,  z) 

=  -^C^x^  -  C4xy  -  iOyy''  +  Wi"(y,  z)  ,  (8) 

I  where  Cp  Bj,  and  are  defined  by  equations  (5)  and  (7).  The  strains  are  found  by 

substituting  the  displacements  (8)  into  the  strain  relations  (2).  The  stresses  then  follow 
directly  using  the  constitutive  relation  (1). 

It  is  of  interest  to  examine  the  strain  which  is 

=  C2Z  +  C3  .  (9) 

Should  the  laminate  be  a  balanced  composite,  i.e.,  the  laminae  are  symmetrically  stacked, 
according  to  composition  and  orientation  with  respect  to  the  midplane  z  =  0,  then  By  = 

0  and  from  equations  (5)  C3  -  0,  which  results  in  a  case  of  pure  bending.  For  the 
opposite  case,  an  unbalanced  composite  exhibits  an  extensional  strain,  C3,  in  bending. 
Such  coupling  effects  are  common  to  laminated  composites. 
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where  the  stresses  exhibit  no  x-dependence,  which  conforms  to  an  earlier  assumption. 
Substituting  for  the  stresses  in  terms  of  displacements  yields  the  field  equations  for  each 
lamina 


„niTTm  ,  „mT,m  ,  .myrn  ,  „m  ym  ^  ^  m  +  =  0 

^  ^  ^26''>VV  ^  ^4S''>ZZ  ^  '^45.’”'>yz  O 


mijm  ,  ^mTjm  .  „m  ^m  .  „ni  ^ni  .  (  in  ,  in  Nyyin  _  n 
C26U,yy  +  ^4511,22  +  ^2  2^  ^  ’ZZ  ^‘^2  3  +  C44JW,y2  ^ 


j.  .,iT>  Miin  j-  +  cPlw”' 

(C36  +  C4s)U,„„  +  (C23  ^  C44;v,y2  4-  C44W,yy  C3  3  W  ,2j 


=  -(c^iCj  +  c?'3Dy  +  2cf,C^) 


The  boundary  conditions  on  the  free  surfaces  normal  to  y  are 


a  m  2™  =  r™  =  0 

"y  'xy  'yz 


and  on  the  free  surfaces  normal  to  z  are 


m  =  2*”  =  m  ^0 
"z  xz  'yz 


For  continuity  at  the  interfaces,  the  boundary  conditions  are 


(u*",vm  wi")  =  w™+l) 


’  'xz’ 


m+1  ^m+l  ^+1) 
’  ^xz  ’  yz  ^ 


respectively. 


It  is  noted  that  the  corner  conditions  are  ambiguous  in  that  there  are  five  possible 
conditions  out  of  which  only  three  can  be  employed  at  any  one  time.  The  remaining  two 
may  or  may  not  be  satisfied  by  the  solution.  Thus,  combinations  may  be  tried  until  some 
satisfying  results  are  achieved. 
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FINITE-DIFFERENCE  SIMULATION 


Function  Representation 

The  mathematical  basis  for  the  finite-difference  method  is  Taylor’s  Series. 
Referring  to  Figure  2,  the  Taylor  Series  expansion  for  a  function  f  at  some  point  y,  z 
about  the  point  (or  node)  I,  J  is 


f(y,  z)  =  f(l,  J)  -H  yf,y(I,  J)  +  zf,^a,  J) 

+  iy^Tyyd,  J)  +  |z2f,2z(I,J)  +  yzf,yz(I,J)  +  ...  .  (15) 

Thus,  for  the  specific  node  I-l,  J,  the  expansion  is 

f(I-l,J)  =  f(l,j)  -  hlf,y  -h  |hiH,yy  '  ...  .  (16) 

Writing  similar  expansions  for  the  remaining  seven  points'  neighboring  the  node  1,  J  and 
simultaneously  solving  expansions  for  the  first  and  second  derivatives  yields  the 
finite-difference  approximations  for  these  derivatives.  All  but  the  last  of  these 
expressions,  given  below,  are  taken  from  Forsythe  and  Wasow  [13] .  They  are 


+  o(h^) 

y  hi  +  hj  Lhj  hi  J  hih2 


J  -1-  1)  -  f(I,  J  -  1) 


f,z(I,  J)  =  ^  J  +  1)  -  fO.  J  -  1)  J  +  0(h") 

f,  (1,  J)  =  -4ir  f(I  +  1,  J)  +  ^  f(I  -  1,  J)]  -  ^  f(I,  J)  +  0(h2) 
yy  hi-l-h2Lh2  hi  J  hihj 

f’zz«’J)  =  +  +  faJ-1)  -  2f(I,J)]  +  0(h2) 


yy  hi-l-h2Lh2  hi 


f>vz(l>  J)  =  ^1,  J  ,  ffd  +  1>  J  +  1)  -  f(I  -  1,  J  +  1)  -  f(I  +  1,  J  -  1) 
y^  2h3(hi  -t-h2)  L 

-I-  f(l-  1,  J  -  1)]  +  0(h2) 
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where  h  is  an  order  of  magnitude  equal  to  hi,  h2,  or  h3.  The  difference  equations  (17) 
are  “central”  differences. 

At  boundaries  and  interfaces  it  is  convenient  to  use  “forward”  and  “backward” 
differences.  Such  difference  equations  are  one-sided  in  that  they  express  a  boundary 
point  in  terms  of  neighboring  points  interior  to  the  boundary.  For  the  present  problem, 
only  first  derivatives  are  of  concern. 

To  derive  such  difference  equations,  expand  two  points,  both  lying  on  one  side  of 
the  reference  point  I,  J,  by  using  equation  (15)  in  conjunction  with  Figure  2.  For 
example,  a  forward  expansion  yields 


f(l+l,J)  =  f(I,  j)  +  h2f,ya  J)  +  J)  +  OCh^^) 

f(I  +  2,  J)  =  f(I,  J)  +  2h2f,y(I,J)  +  i(4h2^)f,yy(I,  J)  +  OCh^")  .  (18) 

Subtracting  one  expression  from  the  other  to  eliminate  the  second  derivative  leads  to  the 
difference  equation  for  the  first  derivative.  Thus,  the  forward  differences  are 

f,y(I,  J)  =  ^  [4f(I+l,J)  -  3f(I,  J)  -  f(I  +  2,J)]  -  0(h2^) 

J)  =  ^  [4f(I,  J  +  1)  -  3f(I,  J)  -  f(I,  J  +  2)]  -  Odia^)  .  (19) 

Similarly,  the  backward  differences  are 

f,y(I,  J)  =  ^  [3f(I,  J)  +  f(I  -  2,  J)  -  4f(I  -  1,  J)]  +  0(hi2) 

J)  =  ^  [3f(I,  J)  +  f(I,  J  -  2)  -  4f(I,  J  -  1)]  +  Odia^)  .  (20) 

It  should  be  pointed  out  that  more  simplified,  but  less  accurate,  forward  and  backward 
expressions  can  be  written;  however,  the  present  application  requires  all  the  accuracy  that 
it  is  possible  to  attain  near  the  free  boundaries.  Thus,  the  higher  order  difference  was 
chosen.  In  addition,  this  choice  yields  a  magnitude  of  error  equal  to  that  found  in 
equations  (17). 
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Using  the  representations  just  obtained,  equations  (11)  through  (14)  can  be 
transformed  into  difference  equations  characterizing  the  problem.  For  example,  the  last 
equation  in  ( 1 1 )  becomes 

(  (ofe  +  [U(I  +  1,  J  +  1)  -  U(I  -  1,  J  +  1)  -  U(I  +  1,  J  -  1) 

2h3(hi  +  hi)  ( 

+  U(I  -  1,  J  -  1)]  +  (cfa  +  0^4)  [V(I  +  1,  J  +  1) 

-  V(I  -  1,  J  +  1)  -  V(I  +  1,  J  -  1)  +  V(I  -  1,  J  -  1)]  j 

+  rW(I  +  1,  J)  +  5^  W(I  -  1,  J)1 

hj  “I"  ]i2  L  hi  -J 

+  ^  cfa  [W(I,  J  +  1)  +  W(I,  J  -  1)] 
ha 

-  2(cf4  +  ^^^cf3)W(I,  J)  =  -hi  hi  [cfiCi  +  cfiDy 

ha 

+  ,  (21) 

where  the  layer  number,  m,  is  left  off  U,  V,  and  W  since  their  location  is  determined  by 
the  node  I,  J. 

Developing  the  Matrix  Equation 

In  this  section,  the  difference  equations,  like  (21),  are  transformed  into  a  linear 
matrix  equation  of  the  form 

[A]  [X]  =  [B]  ,  (22) 

where  A  is  an  N  X  N  coefficient  matrix  (N  being  the  number  of  unknowns  or  equations), 
X  is  the  solution  vector,  and  B  is  the  load  or  input  vector.  To  accomplish  this,  the  three 
unknowns  (U,  V,  and  W)  must  he  uniquely  collapsed  into  the  single  unknown  X  so  that 
at  each  node  three  unique  equations  in  X  will  be  created.  For  instance,  let 
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U^X(l) 

V  ^  X(2)  • 

W->X(3)  ^ 


at  Node  1 


U  ^  X(4) 

V  ^  X(5)  I  at  Node  2  .  (23) 

W  ^  X(6)  ' 


It  remains  to  generalize  such  a  transformation  for  all  nodes. 

It  is  convenient  to  follow  Pipes  [  1  ]  and  his  notation  is  adopted.  If  LAT  is  the 
number  of  nodes  in  one  column  along  the  vertical  axis  (LAminate  Thickness  direction), 
then  the  nodes,  unknowns,  and  equations  can  be  identified  by  a  unique  number  in  terms 
of  the  nodal  position  (I,  J).  If 


JJl  =  3[LAT(I  -  1)  +  J]  -  2 


(24) 


then 


NODE  =  LAT(  I  -  1)  +  J 
U(I,  J)  =  X(JJl) 

(25) 

V(I,  J)  =  X(JJ1  +  1) 

W(I,  J)  =  X(JJ1  +  2) 


Number  the  1st  equation:  JJl 
Number  the  2nd  equation:  JJl  +  1 

Number  the  3rd  equation:  JJl  +  2  .  (26) 


Letting  I  =  1  and  J  =  1 ,  2  consecutively  generates  the  results  in  (23). 

Since  the  finite-difference  equations  involve  unknowns  at  nodes  neighboring  the 
JJl  node,  it  is  necessary  to  develop  transformation  relations  like  (24)  in  order  to  number 
unknowns  at  these  nodes  as  well.  For  example,  using  I,  J  as  the  reference  node,  a 
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transformation  relation  for  an  unknown  at  the  node  I  -  1 ,  J  +  1  is  found  by  letting  I  ->  I 
-  1  and  J  J  +  1  in  (24)  and  giving  the  result  a  unique  name,  for  example  JJ7.  Thus, 

JJ7  =  3[LAT(I-2)  +  J]  +  1  .  (27) 

Using  Table  1,  which  identifies  all  the  unknowns  at  nodes  neighboring  I,  J,  and  following 
the  above  procedure  yields  the  transformation  relations  that  uniquely  number  each 
unknown.  In  summary,  all  of  these  transformations  are 

JJl  =  3*(LAT*I1  +J)  -  2 
JJ2  =  3*(LAT*I2  +  J)  -  2 
JJ3  =  3*(LAT*I2  +  J)  -  5 
JJ4  =  3*(LAT*I  +  J)  -  2 
JJ5  =  3*(LAT*I  +  J)  +  1 
JJ6  =  3=^(LAT*I1  +J)  +  1 
JJ7  =  3*(LAT*I2  +  J)  +  1 
JJ8  =  3*(LAT*I1  +  J)  -  5 
JJ9  =  3*(LAT*I  +  J)  -  5 
JJIO  =  3*(LAT*I1  +  J)  -  8 
JJl  1  =  3*(LAT*(I  +  1)  +  J)  -  2 
JJl  2  =  3*(LAT*I1  +J)  +  4 

JJ13  =  3*(LAT*(I  -  3)  +  J)  -  2  ,  (28) 


where 


11  =  I  -  1 

12  =  I  -  2 


(29) 
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TABLE  1.  NODE  IDENTIFICATION 


— 

Node 

U 

V 

W 

I,  J 

X(JJl) 

X(JJ1  +  1) 

X(JJ1  +  2) 

I  -  1,  J 

X(JJ2) 

X(JJ2  +  1) 

X(JJ2  +  2) 

I  -  1,  J  -  1 

X(JJ3) 

X(JJ3  +  1) 

X(JJ3  +  2) 

I  +  1,  J 

X(JJ4) 

X(JJ4  +  1) 

X(JJ4  +  2) 

I  +  1,  J  +  1 

X(JJ5) 

X(JJ5  +  1) 

X(JJ5  +  2) 

I,  J  +  1 

X(JJ6) 

X(JJ6  +  1) 

X(JJ6  +  2) 

I  -  1,  J  +  1 

X(JJ7) 

X(JJ7  +  1) 

X(JJ7  +  2) 

I,  J  -  1 

X(JJ8) 

X(JJ8  +  1) 

X(JJ8  +  2) 

I  +  1,  J  -  1 

X(JJ9) 

X(JJ9  +  1) 

X(JJ9  +  2) 

I,  J  -  2 

X(JJIO) 

X(JJ10  +  1) 

X(JJ10  +  2) 

I  +  2,  J 

X(JJll) 

X(JJ11  +  1) 

X(JJ11  +  2) 

I,  J  +  2 

X(JJ12) 

X(JJ12  +  1) 

X(JJ12  +  2) 

I  -  2,  J 

X(JJ13) 

X(JJ13  +  1) 

X(JJ13  +  2) 

Generation  of  the  matrix  equation  (22)  now  remains.  To  do  this,  straightforward 
substitution  for  U,  V,  and  W,  using  Table  1,  into  equations  (II)  through  (14)  yields  the 
desired  results  in  equation  form.  For  example,  equation  (21)  becomes 
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ou  ru\u  ^  +^^^5)  [X(JJ5)  -  X(JJ7)  -  X(JJ9)  +  X(JJ3)] 

2h3(hi  +  ha)  ’ 

+  (c™  +  c?\)  [X(JJ5  +  1)  -  X(JJ7  +  1)  -  X(JJ9  +  1) 

+  X(JJ3  +  i)]  }  +  _2hi_  [X(JJ4  +  2)  +  X(JJ2  +  2)] 

’  hi  +  ha  hi 

+  c?l[X(JJ6  +  2)  +  X(JJ8  +  2)] 
hs^ 

-  2(0^4  c?i3)X(JJl  +2) 

ha"' 

=  -hi  ha  [c^3  Ca  +  cfa +  Icf^  C4 1  •  (30) 

To  assure  non-zero  diagonal  terms  in  the  A-matrix,  an  appropriate  equation  number  for 
(30)  is  JQ2  (in  this  case  there  is  only  one  possibility)  where 

JQ2  =  JJl  +  2  .  (31) 

Now,  from  equation  (30),  the  only  nonzero  elements  for  the  JQ2  row  in  the  A-matrix 
are 

A(JQ2,  JJ5)  =  A(JQ2,  JJ3)  =  C 
A(JQ2,  JJ7)  =  A(JQ2,  JJ9)  =  -C 
A(JQ2,  JJ5  +  1)  =  A(JQ2,  JJ3  +  1)  =  D 
A(JQ2,  JJ7  +  1)  =  A(JQ2,  JJ9  +  1)  =  -D 
A(JQ2,  JJ4  +  2)  =  2hic?i/(hi  +ha) 

A(JQ2,  JJ2  +  2)  =  (ha/hi),  •  2hic5‘4/(hi  +  ha) 

A(JQ2,  JJ6  +  2)  =  A(JQ2,  JJ8  +  2)  =  hihaC^^s/hj^ 

A(JQ2,  JJl  -h  2)  =  -2(c?],  +hihacf3/h3') 


(32) 


where 


C  =  hMcf,  +cfs)l2h,(h^  +h2) 

D  =  hihjCc^  +  c54)/2h3(hi  +112)  •  ^^3) 

Note  that  the  material  constants  c5^4  and  c™  are  non-zero  ensuring  a  non-zero 
diagonal  element  A(JQ2,  JJ 1  +  2).  In  addition  to  this,  the  load  vector  is 

B(JQ2)  =  -hih2[c?\C2  +  +  20^604]  •  (34) 


Of  course,  these  results  only  apply  to  node  numbers  where  the  third  equilibrium 
equation  in  (11)  holds.  The  computer  program  logically  connects  appropriate  equations 
with  each  node.  The  matrix  elements  for  the  remaining  equations  (11)  through  (14)  are 
generated  in  a  similar  fashion. 


The  Mesh  Simulation 

The  continuum  is  to  be  simulated  by  a  number  of  nodal  points  that  form  a 
finite-difference  mesh.  The  mesh  is  distributed  over  a  cross  section  of  the  laminate  as 
shown  in  Figure  3.  The  mesh  is  defined  by  the  following  parameters: 

NLAY :  the  number  of  laminae 

LAT:  the  number  of  nodes  along  one  column  in  the  LAminate 

Thickness  direction 

LAW:  the  number  of  nodes  along  one  row  in  the  LAminate  Width 

direction 

FSWl:  the  first  change  in  nodal  spacing  termed  Fine  Simulation 

Width 

K:  magnification  factor  of  the  fine  simulation  width 

H:  the  fine  simulation  width 
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Given  these  parameters,  the  following  parameters  can  be  determined: 


INF(M):  values  of  J  at  the  upper  INterFace  of  the  mth  layer 

FSW2:  the  second  change  in  nodal  spacing 

KH:  the  coarse  simulation  width  (K  =  1,  2,  3,  ...) 

JQMAX  =  3*LAT*LAW:  the  number  of  unknowns  or  equations 
IBW  =  2*(3*LAT  +1):  the  half  bandwidth 

NBAND  =  2*IBW  +  1 :  the  full  band 

The  bandwidth  of  the  coefficient  matrix  is  found  by  considering  that  the  maximum 
number  of  nodes  involved  in  the  difference  equations  is  three,  as  can  be  seen  from 
expressions  (19)  and  (20),  and  calculating  the  maximum  number  of  consecutive  elements 
on  both  sides  of  the  diagonal  to  and  including  the  last  off-diagonal  non-zero  element. 

Selecting  equations  representing  the  conditions  to  be  imposed  at  each  node 
remains  to  be  accomplished.  Because  of  the  arbitrariness  of  the  corner  conditions,  a 
number  of  choices  are  possible.  Those  selected  for  this  program  are  illustrated  in  Figure 

4. 

A  user’s  guide  and  a  more  detailed  description  of  the  computer  program  are 
presented  in  Appendix  C.  A  program  listing  is  provided  also  in  Appendix  C. 

RESULTS 


The  results  given  below  were  obtained  using  a  square  mesh,  magnification  factor 
K  =  1,  of  size  (LAW,  LAT)  =  (13,  9).  A  complete  mesh  description,  taken  from  the 
program  output,  is  displayed  in  Table  2.  It  is  seen  that  these  dimensions  represent  a  beam 
rather  than  a  plate.  The  program  was  run  on  an  IBM  370  computer  utilizing  virtual 
storage. 

A  single  material  having  properties  typical  of  a  high  modulus  graphite-epoxy  was 
chosen  for  the  above  mesh.  Using  standard  notation. 

Ell  =  20.0  X  10^  psi  ,  ^^12  ^13  -  ^23  =  0-21 

E22  ~  E33  =  2.1  X  10^  psi 


Gi2  ~  Gi3  =  G23  =  0.85  X  10^  psi 
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TABLE  2.  MESH  DESCRIPTION  TAKEN  FROM  PROGRAM  OUTPUT 


8ENTIHG  -OF  A  LAMINATFO  PLATE  *>!'*- 


INPUT  CAFA 

- - - vjM  -i  £  R-  t>F  -  LAYERS  INC  RUS  S . S  FC  T  f  CH?-  NtAY— ^ - A- 

.  NU'^  TER  OF  MODES  ON  VE  RT  I  CAi  AX  f  St  t  AT - =  -13  ' 

. - . -  . -  NL-AMcR  OF  HOOFS  ON  HORfZONTAL-  --AXTS>  1^0-= - P- 

C^5^Gc  IM  '^FSH  WIDTH  (FSWl)  AT  I  =  3 

CHivGF  IN  MESH  WIDTH  (FSW2)  AT  I  =  7 

MESH  WIDTH  MAGNIFICATION  FACTOR,  K  =  I 


L  AVE  I  NO.  I  INTERFACE  AT  J  =  4 


- tr  A  Y— R“  fifO . P — — f-NT  FR-  F  A  C  F  — AT— J- 


LAYET  NO.  3  INTERFACE  AT  J  =  10 

- NO-:^- - 4 - PFFFRFACE  AT  J  =— TO - 

- ALiON  -W-i  OTH,--  -H--  = . 0-.:0<n-6-7- - 


where  the  subscript  “l”  refers  to  the  fiber  direction.  The  two  laminate  configurations 
which  are  considered  are 

A  =  [d,  0,0,0] 

and 

B  =  [0,  0,0,0] 
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with  6  as  in  Figure  1  such  that  0  degree  <  ^  <  90  degrees.  Typical  laminate  data  and 
load  constants  are  displayed  in  Table  3.^  Here  the  additional  constant  MT  is  the  resulting 
moment  required  to  produce  a  specified  maximum  strain  which,  for  the  present  analysis, 
is  =  1.0  X  10"^  inch/inch  (see  Appendix  B). 

A  sample  of  the  results  for  the  displacement  functions  U,  V,  and  W  is  presented 
in  Table  4.  Examination  of  their  variation  with  respect  to  z  reveals  the  apparent 
symmetry  relations. 


U,  V  antisymmetric  in  z 

W  symmetric  in  z 


within  an  accuracy  of  two  digits. 

Symmetries  with  respect  to  y  are  evident  for  the  strains  within  three-digit 
accuracy.  Samples  of  these  results  are  plotted  in  Figures  5  and  6.  Coupling  these  apparent 
symmetries  with  the  strain  relations  i2)  in  an  expanded  form  yields 


U,  V  antisymmetric  in  y 

W  symmetric  in  y 


The  displacement  results  verify  this  precisely  for  U  (to  four  places),  but  show  some 
deviation  in  V  and 

To  illustrate  the  effect  of  bending  on  the  stress  distribution.  Figures  7  through  19 
are  presented.  Although  convergence  to  the  exact  values  has  yet  to  be  demonstrated,  the 
results  do  have  quahtative  merit.  The  following  cases  result  from  a  bending  strain  of  e^  = 
1,0  X  10"^  inch/inch  prescribed  at  the  bottom  surface. 

Of  principal  interest  are  the  interlaminar  stresses  illustrated  in  Figures  7  through 
12.  We  note  that  laminates  composed  of  30  degree  or  45  degree  layers  produce  the 
greatest  stress  rise  in  at  the  free  edge  with  a  more  pronounced  effect  occurring  if  the 

angle  plies  are  on  the  outside,  i.e.,  system  A  =  [6,  0,  0,  0] .  A  similar  effect  is  seen  in  the 
shear  stress  Ty^,  although  the  rise  in  stress  is  sharply  blunted  by  the  requirement  of  zero 

3.  The  thermal  problem  is  neglected  in  this  preliminary  analysis  even  though  expansion 
coefficients  appear  in  the  program. 

4.  It  is  interesting  to  note  that  the  y-symmetries  for  V  and  W  are  verified  precisely  using  the 
coarser  mesh  (LAW,  LAT)  =  (8,  9)  which  decreases  the  relative  size  of  the  bandwidth. 
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TABLE  3.  (Concluded) 


TABLE  4.  DISPLACEMENT  FUNCTION  RESULTS  TAKEN  FROM 
PROGRAM  OUTPUT  FOR  LAMINATE  DESCRIBED  IN 
TABLES  2  AND  3 


. 

"9--FSPL  AC 

NC)E 

)-m  splacemfnt 

V-nrSPLACGMENT 

W-DFSPLACEMENT 

1 

0. 1615610-04 

0.2646360- 34 

-0.9090390-05 

- . 0,  1495  800— 04 

0,2198310-04  - 

-0. 9353040-05 

0.1253810-04 

0.173176D-Q4 

-0. 9519390-05 

. .  v.-\  jric:  "a  c  O  A  c .  - 

- - - A 1-?  7  0 1 4  n  -  04 . - 

-0. 953590D-05 

5 

0.6116060-05 

0. 318957D- 35 

-0, 9400020-05 

. 4.- . 

-  0.  0043550-05 

- . -  O.403364D-O5- 

- 0.924686D-05 

1 

0. 4871390-09 

0. 2507690-08 

-0.9165890-05 

c 

A  rwn^i-  — 

-  0  *40  391  30^05 - 

- 0.9250840-05  - 

Tif .  j 

i 

-0.61 15750-05 

-0.3194323-35 

-0.9376930-05 

. — . l-l . 

- .  -0* 0?66B9D-05  . - 

-0.1263080-04 

-0.9528  920-05 

11 

-0.1253540-04 

-0.1732110-04 

-0.9544130-05 

.  -  -  \  2  . - 

.  -0.  1495  50D-04- . 

-  -0.219  3309-04- 

.  -0.9371610-05 

1  3 

-0. 1615030-04 

-0.2643030- 04 

-0.9098910-05 

stress  at  the  free  edge,  and  here  the  stress  in  system  B  =  [0,  0,  0]  is  slightly  more 

pronounced  than  that  in  A.  The  largest  stress  rise,  an  order  of  magnitude  greater  than 

and  Ty^,  is  created  in  the  A-system  in  Again  it  is  the  30  degree  laminate  incurring 

the  sharpest  stress  rise,  but  here  the  15  degree  laminate  overshadows  the  45  degree 
laminate.  In  summary,  the  laminates  containing  15  degree  through  45  degree  layers 
located  adjacent  to  0  degree  layers  have  the  largest  interlaminar  stresses  for  the  cases 
considered;  i.e.,  0  degree  <  0  <  90  degrees  taken  in  15  degree  intervals. 

Some  results  peculiar  to  the  numerical  method  of  solution  should  be  pointed  out. 
Referring  to  Figure  9,  we  note  a  sharp  rise  in  the  stress  at  the  midpoint  node  (I,  J)  = 

(5,  7).  This  is  a  result  of  fixing  the  displacements  at  I  =  5  and  6,  J  =  7  in  the  program  in 
order  to  zero-out  rigid  body  motion  and  drift  in  the  solution  routine.  However  averaging 
the  values  for  just  above  and  just  below  the  interface  (at  J  =  7,  m  =  2  and  m  =  3) 

yields  a  more  plausible  result.  Since  the  tractions  must  be  continuous  at  the  interface 
anyway,  this  averaging  technique  was  also  applied  at  the  free  edges  where  the  free  surface 
conditions  were  adopted  in  heu  of  the  continuity  conditions.  This  technique  had  varying 
success  as  illustrated  by  the  75  degree  and  90  degree  configurations  in  Figures  10  and  11. 

The  in-plane  stresses  are  illustrated ' in  Figures  13  through  19.  In  Figure  13,  we 
find  that  in  the  0  degree  layers  is  independent  of  the  orientation  of  the  adjacent  layer 

when  the  maximum  strain  is  specified.^  This  facilitates  the  presentation  of  both  systems 
A  and  B  in  one  figure.  It  is  interesting  to  note  in  Figure  14  that  rises  at  the  free  edge 

if  the  0  degree  layers  are  outside  the  laminate  and  drops  if  these  layers  are  inside  the 
laminate. 

Observation  of  Figures  1 5  and  1 7  for  the  distribution  of  Uy  and  r^^y  with  respect 

to  z  reveals  that  the  off-axis  layers,  particularly  again  for  15  degrees  through  45  degrees, 
serve  as  stress  raisers  with  the  effect  considerably  more  pronounced  if  the  0  degree  layers 
are  inside. 

Typical  distributions  of  Uy  and  r^y  with  respect  to  y  are  shown  in  Figures  18  and 

19.  The  disturbing  feature  of  these  plots  is  that  the  stresses  just  above  an  interface  do 
not  approach  zero  at  the  free  surface.  One  cause  of  this  problem  is  the  placement  of 
nodes  directly  on  the  interface,  which  requires  their  occupation  by  both  layers.  Then  at 
the  corners,  as  stated  previously,  the  multitude  of  boundary  conditions  cannot  be 
satisfied.^  However  this  problem  is  confined  to  the  free  surface  nodes  and  one  line  of 


5.  In  agreement  with  the  beam  theory  approximation. 

6.  Placing  the  interface  between  two  nodal  lines  may  alleviate  this  problem. 
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interior  nodes.  To  see  this,  one  may  examine  the  curves  for  the  A-system  at  J  =  4-  and  J 
=  10+  and  note  that  they  are  reflections  of  each  other  within  the  range  3  <  I  <  7.  Since, 
from  above,  Oy  and  appear,  in  general,  to  be  antisymmetric  in  z,  the  correct  values  at 

J  =  10+  are  recovered  within  this  range  if  we  accept  the  values  at  J  =  4-. 


CONCLUSIONS 


Although  only  two  types  of  laminate  systems  were  considered,  namely  A  =  [0,  0, 
0,  6]  and  B  =  [0,  6,  6,  0],  it  is  reasonable  to  assume  from  these  results  and  from 
physical  considerations  that  the  following  symmetry  relations  hold  for  balanced  (Bjj  =  0) 
composites: 


U,  V  antisymmetric  in  y  and  z 
W  symmetric  in  y  and  z 


where  U,  V,  and  W  are  displacement  functions  of  y  and  z.  Based  on  the  stress  results, 
laminates  containing  layers  oriented  within  the  range  15  degrees  <  0  <  45  degrees 
produce  the  largest  interlaminar  stresses  out  of  the  cases  studied,  0  degrees  <  0  <  90 
degrees  taken  in  1 5  degree  intervals.  In  fact  this  same  group  of  laminates  produces  high 
values  in  the  in-plane  stresses  as  well,  with  the  effect  considerably  more  pronounced  for 
the  A-system.  Although  some  deviations  in  stress  occur  in  the  numerical  solution,  they 
are  localized  to  a  double  line  of  nodes  at  the  boundary.  This  is  a  disconcerting  feature  of 
the  solution  in  that  the  boundary  region  stresses  appear  to  be  critically  involved  in 
delamination-type  failure,  which  makes  their  accurate  determination  desirable. 

This  study  provides  a  base  for  future  work  in  this  area.  Using  the  present  program 
coupled  with  an  out-of-core  equation  solver  routine,  unbalanced  laminates  may  be 
studied.  Using  the  symmetry  relations  discussed  above,  the  present  computer  program 
may  be  modified  to  more  efficiently  handle  balanced  laminates  (By  =  0). 
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Figure  2.  Finite-difference  mesh. 


24 


Figure  4.  Equations  selected  for  each  node. 
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Figure  14.  Variation  of  the  normal  stress  (symmetric  in  y)  with  y 
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Figure  17.  Variation  of  the  shear  stress  with  z. 


APPENDIX  A 


LAMINATE  CONSTANTS 


Following  Reference  9  or  10,  define 


^33 


m 


i,  j  =  1,  2,  6 


and  let  t  be  the  half-thickness  of  the  laminate,  ho  the  thickness  of  a  lamina,  and  N  the 
total  number  of  layers;  then 


N 

nj  -^0 

m=l 


A.-  =  h, 


m 


Vi  2 


N  N 

Qij"'(2m-l)  -  N  Qi^ 

m=l  m=l 


h  ^ 

D'-i  =  ^ 

«  3 


r  N 

Qjj™(3m^  -  3m  +  1) 

I  m=l 


N 

|n  I]  Qij”’(2m-1) 

m=l 


N 


+  I  N=  Yi  Qii“ 

m=l 


with  i,  j  =  1 ,  2,  6.  Finally,  let 


A*  =  A'*  ,  B*  =  -A'^B  ,  and  D*  =  D  -  BA‘‘B 
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where  the  letters  symbohze  3X3  matrices.  Then, 
B'  =  B*(D*)'^ 


and 


D'  =  (D*r 


Considerable  simplification  is  attained  if  the  laminate  is  balanced,  which  implies  Bjj  =  B 

=  0. 


ij 


APPENDIX  B 


STRAIN  SPECIFICATION 


Rather  than  prescribe  the  laminate"  loading  as  end  moments,  the  maximum  strain, 
e  at  the  top  and  bottom  surfaces,  z  =  will  be  prescribed.  From  equation  (9), 

we  have 

max  =  c^z^iax  + 

A. 

Now  from  equations  (5), 


C3  =  -BiiM  = 

A-'l  1 


SO  that 


max  ^ 


=  C. 


,max 


D'li 


and,  thus. 


C2 


D'li 


B'l  1  +  D'l  1 


max 


^max 


In  the  computer  program,  we  set  -  -1.0  X  10  ^  inch/inch  at  the  top  surface  z 

-(.zmax  evaluate  the  constant  Cj  which  represents  the  inverse  bending  radius. 
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APPENDIX  C 


THE  COMPUTER  PROGRAM 


Program  Description 

The  computer  program  is  an  in-core  program  and  is  not  overlayed.  It  is  felt  that  a 
flow  chart  of  the  program  would  be  no  less  complicated  than  the  presentation  of  a  listing 
with  an  accompanying  explanation,  so  the  latter  choice  will  be  followed.  Certain 
statements  in  the  program  are  extraneous  to  the  problem  in  this  report  because  the 
program  is  in  steady  transition  to  handle  more  general  problems.  A  part-by-part 
description  follows. 

Part  L  Part  I  contains  a  brief  definition  of  terms  and  an  explanation  of  the  order 
and  format  of  the  data  cards.  The  dimensions  of  the  data  are:  H  is  in  inches,  E  is 
material  constants  in  psi  (the  shear  moduli  G12,  etc.,  are  read  into  the  El 2,  etc.,  arrays), 
ALPHA  is  the  coefficient  of  expansion  in  inches/inch/^^F,  and  THETA  is  the  lamina 
orientation  in  angular  degrees.  Precision  and  dimension  statements  are  then  established, 
data  are  read  in,  and  mesh  parameters  are'  calculated.  The  letter  M  refers  to  the  layer 
number.  In  the  loop,  DO  9000,  IRAN  counts  each  laminate  layup  from  one  to  IRUN 
(only  changes  in  lamina  orientation  are  allowed  for  within  this  loop). 

Part  II.  Part  II  calculates  the  anisotropic  stiffness  matrix.  BETA  is  in  radians. 
CPU,  etc.,  are  the  orthotropic  elastic  constants  in  the  primed  coordinate  system. 
Cll(M),  etc.,  are  the  anisotropic  elastic  constants  for  the  Mth  lamina  in  the  x,  y,  z 
coordinate  system.  ALIP(M),  etc.,  are  the  coefficients  of  thermal  expansion  in  the 
primed  coordinate  system  and  ALl(M)  are  those  coefficients  in  the  x,  y,  z  coordinate 
system,  both  the  the  Mth  lamina.  Finally,  the  subroutine  MATCON,  which  calculates  the 
laminate  MATerial  CONstants,  is  called. 

Part  III.  Part  III  calculates  the  coefficient  matrix  for  the  difference  equations. 
The  loops  DO  100  and  DO  101  count  through  the  mesh  node-by-node.  DO  3000  zeroes 
out  the  A-matrix. 

The  logic  that  associates  the  various  field  conditions  with  each  node  and  correctly 
fills  out  the  A-matrix  is  contained  in  DO  102.  First  the  node  I,  J  is  tested  to  determine 
the  proper  layer  number,  M.  Then  the  node  is  checked  to  see  if  it  lies  on  a  boundary, 
along  J  equals  a  constant,  or  lies  at  some  select  position  (in  this  case,  IMID  or  JMID).  If 
it  does,  the  program  is  routed  to  the  statement  number  that  contains  the  non-zero  matrix 
elements  satisfying  the  conditions  imposed  at  this  node.  Should  the  node  not  lie  at  any 
of  these  preselected  locations,  the  program  passes  through  the  IF  statements  on  J  to 
statement  number  193,  which  initiates  a  series  of  checks  to  see  if  the  node  lines  on 
selected  values  of  I.  These  values  include  the  boundaries  I  =  1  or  I  =  LAW,  the  changes  in 
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nodal  spacing  I  =  FSWl  or  I  =  FSW2,  and  all  points  in  the  region  between  FSWl  and 
FSW2.  Should  the  node  not  lie  at  any  of  these  locations,  the  program  passes  through  the 
IF  statements  on  I  and  evaluates  the  non-zero  coefficients  for  the  only  remaining 
possibility,  the  equilibrium  terms  for  a  square  mesh. 

When  a  node  does  lie  on  some  select  location,  say  J  equals  LAT,  then  the  logic  in 
that  statement  series,  say  the  series  starting  from  statement  number  202,  guides  the 
program  through  the  checks  on  selected  values  of  I  in  a  fashion  similar  to  that  above. 
The  logic  is  easily  understood  by  reading  directly  from  the  listing. 

Upon  reaching  statement  number  102,  the  A-matrix  (3  X  JQMAX)  is  full.  The 
elements  of  the  A-matrix  lying  within  the  bandwidth  are  then  stored  in  the  banded 
matrix  AX.  The  loops  DO  100  or  DO  101  then  continue  for  the  next  node,  if  any.  The 
previous  A-matrix  is  destroyed  and  regenerated  for  the  new  node  until  the  loop  DO  100 
is  satisfied. 

At  rewind  9,  the  matrix  AX  and  the  load  vector  X  are  stored  for  later  use.  The 
loop  DO  107  stores  the  load  vector  X(I)  in  AX(I,  NBD).  Then  a  series  of  WRITE 
statements  (listed  as  comments)  will  output  the  coefficient  matrix  AX  and  load  vector  X 
should  they  be  desired.  Finally,  the  solver  routine,  TRMSTR,  is  called.'^ 

Part  IV.  Part  IV  outputs  the  functional  displacements  and  provides  an  accuracy 
check.  Just  below  statement  number  4006,  the  STOP  1  statement  will  terminate  the 
program  if  the  coefficient  matrix  AX  is  singular.  (Such  an  occurrence  probably  indicates 
an  error.)  The  loop  DO  108  stores  the  solution  vector  AX(I,  1)  in  X(I).  Then  the  original 
values  of  the  matrix  AX  and  load  vector  X  are  read  back  into  the  AX  array  and  R 
vector,  respectively. 

The  loops  DO  11  and  DO  12  output  the  values  for  the  functions  U(y,  z),  V(y,  z), 
and  W(y,  z)  which  occur  in  the  displacements  u,  v,  and  w,  respectively. 

The  series  of  statements  from  the  one  above  9950  to  9990  outputs  the  accuracy 
results.  These  results  provide  the  difference  between  the  original  load  vector,  now  stored 
in  the  R-array,  with  the  calculated  load  vector,  which  is  found  by  substituting  the 
appropriate  solution  vectors,  X(I),  into  each  matrix  equation.  In  addition  to  giving  the 
accuracy  of  each  equation,  an  average  accumulated  accuracy  is  provided. 

Part  V.  Part  V  outputs  the  strains  and  stresses.  The  logic  is  similar  to  that  in  Part 
III.  Knowing  that  the  finite-difference  relations  for  the  strains  differ  for  various  mesh 
locations,  the  strains  are  split  into  terms  dependent  upon  the  value  of  I  and  terms 
dependent  upon  the  value  of  J.  The  strain  SX,  which  represents  e.^,  depends  upon  neither 

the  value  of  I  nor  the  value  of  J  and  is  determined  prior  to  any  logical  branching. 


7.  Actually  the  AX-matrix  stores  a  transposed  A-matrix;  i.e.,  instead  of  storing  row  elements 
crosswise  or  in  a  row,  they  are  stored  in  the  AX-matrix  vertically  or  in  a  column.  The  result 
is  a  drastic  reduction  in  “wall-time”  on  the  IBM  370.  This  necessitated  a  slight  revision  in 
the  solver  routine,  TRIMSS,  as  written  by  Billy  Gibbs,  U.S.  Army,  Redstone  Arsenal  [  14] . 
So  here  it  is  called  TRMSTR  or  TRIMSS  transposed. 
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First,  the  node  is  checked  to  determine  its  location  with  respect  to  I,  and 
I-dependent  strains  (or  the  partial  strain,  SYZI)  are  calculated.  Then  the  loop  DO  392 
establishes  the  correct  layer  number,  M,  in  order  to  check  if  J  lies  on  the  interface, 
INF(M).  Upon  determining  the  correct  location  of  the  node  with  respect  to  J,  the 
J-dependent  strains  (or  the  partial  strain,  SYZJ)  are  calculated.  Statement  number  391 
totals  the  partial  strains  to  obtain  SYZ.  The  stresses  are  then  calculated  in  a  straight 
forward  manner  using  equation  (1).  Note  that  the  stresses  are  calculated  twice  at 
interface  nodes,  once  for  the  material  below  the  interface  and  again  for  the  material 
above. 


Part  VI.  The  subroutine  MATCON  calculates  the  MATerial  CONstants  C.-,  BU, 
BV,  and  DV  as  defined  earlier  in  the  text. 

Part  VII.  The  subroutine  MAMULT  is  a  MAtrix  MULTipUer  and  is  easily 
understood  from  the  listing. 

Part  VIII.  The  subroutine  MATIN4  is  a  MATrix  INversion  routine  which  is 
described  in  Reference  14. 

Part  IX.  The  subroutine  TRMSTR  is  the  equation  solver  which  is  described  in 
the  Usting. 

Part  X.  The  subroutine  RITE  is  used  to  wRITE  out  a  matrix  or  vector. 


Program  Listing 

The  complete  listing  of  the  program  is  contained  in  the  following  pages. 
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FORTRAN  IV  G1  RELEASE  2.0 


MAIN 


DATE  =  75007 


08/16/07 


0001 

0002 

0003 

0004 

0005 


0006 

0007 


0008 


0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 


JOMAX  IS  THE  NUMBER  OF  UNKNOWNS  OR  EQUATIONS  TO  BE  SOLVED. 

A  IS  A  FULL  MATRIX  (3  X  JQMAX)  REPRESENTING  EACH  NODE. 

AX  IS  THE  BANDED  MATRIX  <NBAND+1  X  JQMAX). 

X  IS  THE  LOAD  VECTOR.  AFTER  TRIMSS  X  BECOMES  THE  SOLUTION  VECTOR. 

IF  THE  NUMBER  OF  LAYERS  EXCEED  6»  THE  COMMON  /MC/  AND  DIMENSION  (Ell 
E22,  ETC.)  STATEMENTS  MUST  BE  REDIMENSIONED  TO  AGREE  WITH  LAT. 
REMEMBER  TO  PLACE  A  COMMON  /MC/  STATEMENT  IN  SUBROUTINE  MATCON. 


USE  THE  FOLLOWING  ORDER  FOR  DATA  CARDS 


DATA  CARO  NO. 
1 


DATA 

NLAY,  LAT,  LAW,  FSWl,  K 
H 

Ell,  E22,  E33,  £12,  E13,  E23 
NU12,  NU13,  NU23 

ALPHA  1  PRIME,  ALPHA  2  PRIME,  ALPHA  3  PRIME 


FORMAT 

5110 

G12.5 

8G12.5 

8G12.5 

8G12.5 


NOTE,  REPEAT  CAROS  OF  THE  TYPE  3,  4,  5  FOR  EACH  ADDITIONAL  LAYER 

6  SXMAX,  C3E  10G10.3 

7  IRUN  5110 

THETA(l),  THETA(2),  THETA(3),  ETC.  10G10.3 


NOTE,  REPEAT  CARD  8  FOR  EACH  ADDITIONAL  LAYUP. 


INTEGER  P,  FSWl,  FSW2 

DOUBLE  PRECISION  TEST,  R,  ERR,  AVE,  OT 
DOUBLE  PRECISION  AX,  X 
DOUBLE  PRECISION  THETA,  BETA 

DOUBLE  PRECISION  CM,  CN,  CM4,  CN4,  CM3N,  CN3M,  CM2,  CN2,  GNU21, 

1  GNU31,  GNU32,  DET,  CPU,  CP22,  CP33,  CP12,  CP13 

2  CP23,  CP44,  CP55,  CP66 

DIMENSION  AX(162,351),A(3,351),  X{35U,  R(351) 

COMMON  /MC/  Cll(6) ,C12(6) ,C16( 6),C22(6) ,C26(6),C66(6) ,C13(6) , 

1  C23(6) ,C36(6) ,C44(6) ,C45(6) ,C55 ( 6 ) ,C33 ( 6 ) , ALl ( 6 ) , AL2 ( 6 ) , 

2  AL3( 6) ,AL6( 6} ,C2,C3,C3E,C4,BU,DU,BV,DV,H,SXMAX,NLAY, INF (6) 

DIMENSION  Eli( 6) , E22( 6) , E33I 6) ,E12(6),E13{6), E23 ( 6 ) , GNU  1 2 ( 6 ) , 

1  GNU13(6) ,GNU23(6),THETA{6),  AL1P(6),  AL2P(6),  AL3P(6) 


WRITE(6,600) 

READ{5,601)NLAY,LAT,LAW,FSW1,K 

FSW2=LAW-FSW1+1 
JQMAX  =  3*LAW=t=LAT 
IBW  =  2^(3*LAT+1) 

IBWl  =  IBW+1 
NBAND  =  2*IBW+1 

WRITE! 6, 602) NLAY, LAT, LAW, FSWl, FSW2,K 

LATl=LAT-l 

IMID  =  (LAW+l)/2 

JMID  =  (LAT+l)/2 

DO  501  M=l,  NLAY 


00000010 

00000020 

00000030 

00000040 

00000050 

♦00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250' 

00000260 

00000270 

00000280 

00000290 

,00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000580 
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0022 

INF(M)  =  l+M’i‘LATl/NLAY 

00000590 

0023 

WRITE{6,608)M, INF(M) 

00000600 

0024 

501  CONTINUE 

00000610 

C 

00000620 

c 

NOTE  THAT  INF{NLAY)  EQUALS  LAT  AND  IS  NOT  AN  ACTUAL  INTERFACE. 

00000630 

c 

00000640 

0025 

READ(5,603)  H 

00000650 

0026 

WRITE(6t607)  H 

00000660 

0027 

HSQ  =  H**2 

00000670 

c 

00000680 

0028 

WRITE{6,604) 

00000690 

c 

00000700 

0029 

DO  500  M=1,NLAY 

00000710 

0030 

READ(5,603)E11(M) ,E22(M) tE33(M),E12(M)  ,E13(M)  , 

E23(M) 

00000720 

0031 

READ( 5,603 }GNU12(M) ,GNU13(M) ,GNU231M) 

00000730 

0032 

WRITE(6,605)  M,  Ell(M),  E22(M),  E33(M),  E12(M) 

,  E13(M),  E23(M), 

00000740 

1  GNU12(M),  GNU13(M),  GNU23(M) 

00000750 

0033 

READ(5,603) AL1P{M) ,  AL2P(M),  AL3P(M) 

00000760 

0034 

500  CONTINUE 

00000770 

c 

00000780 

0035 

READ! 5,606)  SXMAX,  C3E 

00000790 

0036 

READ{5,601)  IRUN 

00000800 

c 

00000810 

0037 

DO  9000  IRAN  =  1,  IRUN 

00000820 

0038 

READ(5,606)  ( THETA ( M ), M= 1 , NLA Y ) 

00000830 

c 

00000840 

C 

00000860 

C 

CALCULATION  OF  ANISOTROPIC  STIFFNESS  MATRIX  TERMS 

REFERRED  TO  X,Y,Z 

00000870 

C 

00000880 

C 

00000900 

0039 

WRITE(6,613) 

00000910 

0040 

XX  =  0.0 

00000920 

0041 

00  3001  M=1,NLAY 

00000930 

0042 

BETA  =  .0174532925199433D0*THETA(M) 

00000940 

0043 

CM=DCOS(BETA) 

00000950 

0044 

CN=OSIN(BETA) 

00000960 

0045 

IFIDABSICM) .LT.l.E-08)  CM  =  0. 

00000970 

0046 

IF(DABS(CN).LT.l-E-08)  CN  =  0. 

00000980 

0047 

CM4=CM##4 

00000990 

0048 

CN4=CN^4=4 

00001000 

0049 

CM3N=CM««3*CN 

00001010 

0050 

CN3M=CN««3*CM 

00001020 

0051 

CM2-CM*#2 

00001030 

0052 

CN2=CN**2 

00001040 

0053 

GNU21=GNU12{M)*E22(M) /Ell (M) 

00001050 

0054 

GNU31  =  GNU13(M)«E33(  M)  /  ElKM) 

00001060 

0055 

GNU32=GNU23(M)*E33(M) /E22 (M) 

00001070 

0056 

DeT=l.-GNU12(M)>S‘GNU21-GNU23(M)»!'GNU32-GNU13(M}*GNU31 

00001080 

1-2.*GNU12(M)»!'GNU23(M}«GNU31 

00001090 

0057 

CPll^ElKMI’tcJ  l.-GNU23(M)=«'GNU32)/DET 

00001100 

0058 

CP22  =  E22(M)=}‘(  l.-GNU13(M)’!'GNU31)/DET 

00001110 

0059 

CP33=E33(M)*( 1,^GNU12(M)*GNU21 )/DET 

00001120 

0060 

CP12  =  E11(M)«(  GNU21+GNU23(M)’i'GNU31)/DET 

00001130 

0061 

CP13=E11(M)*(GNU31+GNU21*GNU32)/DET 

00001140 

0062 

CP23=E22(M)=<'(GNU32  +  GNU12(M)=*GNU31)/DET 

00001150 

0063 

CP44=E23{M) 

00001160 

49 


FORTRAN  IV  G1  RELEASE  2.0 


MAIN 


DATE  =  75007 


08/16/07 


0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 


0079 

0080 

0081 

0082 

0083 


0084 

0085 

0086 

0087 

0088 

0089 


0090 

0091 

0092 

0093 

0094 


CP55=E13(M)  00001170 

CP66=E12(M)  00001180 

C11(M)=CM4*CP11  +  2.4CM2’S'CN2’«'CP12+CN4*CP22+4.*CM2«CN2*CP66  00001190 

C12(M)=CM2*CN2’»=CP11+(CM4+CN4)*CP12+CM2=<'CN2*CP22-'CM2*CN2=S‘4.*CP66  OOOOiaOO 
C16{M)  =CM3N’!'CPll-(  CM3N-CN3M)=«=CP12-CN3M*CP22-2.*(CM3N~CN3M)><'CP66  00001210 

C22{  M)=CN4>>=CP11  +  2.=5'CM2>!‘CN2«CP12  +  CM4’«'CP22  +  4.«CM2’«‘CN2*CP66  00001220 

C26(M)  =CN3M=i‘CPl  I- ( CN3M-CM3N )  «C PI 2-CM3N’4'C P22-2 <  CN3M'CM3N  ) ’i'C P66  00001230 

C66(M)=CM2*CN2’^CPll‘‘2.#CM2*CN24CP12+CM2=4'CN2*CP22+(CM2“CN2)«=i'2*CP6600001240 
C13{M)=CM2*CP13+CN2*CP23  00001250 

C23(M)=CN2>«^CP13+CM2*CP23  00001260 

C36{M) =CM*CN*(CP13-CP23)  00001270 

C44( M)=CM2#CP44+CN2*CP55  00001280 

C45{M) =CM*CN^(CP55-CP44)  00001290 

C55(M)=CN2*CP44+CM2«CP55  00001300 

C33(M)=CP33  00001310 

r  00001320 

Q  00001330 

(.  00001350 

C  CALCULATION  OF  THE  COEF.  OF  THERMAL  EXPANSION  REFERRED  TO  X,Y»Z  00001360 

Q  00001370 

00001390 

AL1(M)=CM2’>^AL1P(M)+CN2>{'AL2P(M)  00001400 

AL2{M)=CN2*AL1P(M)+CM2*AL2P< M)  00001410 

AL3(M)=AL3P(M)  00001420 

AL6(M)=2.’8‘CM«CN’S‘{AL1P(M)-AL2P(M)  )  00001430 


WRITEi6,620)  M,  Cll(M),  C12(M)t  C13tM),  XX,  XX,  C16(M),  CPU, 

1  CPU,  CP13,  XX,  XX,  XX,  C22(M),  C23(M),  XX,  XX, 

2  C26(M),  CP22,  CP23,  XX,  XX,  XX,  C33(M),  XX,  XX, 

3  C36{M),  CP33,  XX,  XX,  XX,  THETA(M),  C44(M),  C45(M), 


00001390 

00001400 

00001410 

00001420 

00001430 

00001440 

00001450 

00001460 

00001470 

00001480 


C 

C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


WRITE(6,611) 


4  XX,  CP44,  XX,  XX,  C55(M),  XX,  CP55,  XX,  C66<M),  CP66  00001490 

00001500 

3001  CONTINUE  00001510 

00001520 

WRITE(6,611)  00001530 

00001540 

DO  503  M=1,NLAY  00001550 

WRITE{6,614)  M,  THETA(M),  ALUM),  AL2(M),  AL3{M),  AL6(M),  00001560 

1  ALIP(M),  AL2PIM),  AL3P(M)  00001570 

503  CONTINUE  00001580 

00001590 

CALL  MATCON  00001600 

00001610 

00001620 

00001640 

CALCULATION  OF  THE  COEFFICIENT  MATRIX  FOR  THE  DIFFERENCE  EQUATIONS  00001650 

00001660 

00001680 


DO  503  M=1,NLAY 

WRITE{6,614)  M,  THETA(M),  ALUM),  AL2(M),  AL3{M),  AL6(M), 
1  ALIP(M),  AL2PIM),  AL3P(M) 

503  CONTINUE 

CALL  MATCON 


KJl  =  1 

KQl  =  KJl  +  1 

KQ2  =  KJl  +  2 


00001690 

00001700 

00001710 

00001720 


DO  100  1=1, LAW 
DO  101  J=l,  LAT 


00001730 

00001740 
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0095 

C 

00  3000  IM  =  KJl,  KQ2 

00001750 

00001760 

0096 

DO  3000  JM  =  1,  JQMAX 

00001770 

0097 

A(IMtJM)  =  0. 

00001780 

0098 

3000  CONTINUE 

00001790 

0099 

C 

11=1-1 

00001800 

00001810 

0100 

12=1-2 

00001820 

0101 

Z  =  {FLOAT(  J)-(FL0AT(LATH‘1.)/2.)*H 

00001830 

0102 

NODE  =  LAT*I1+J 

00001840 

0103 

JJl  =  3*(LAT*Il+J)-2 

00001850 

0104 

JJ2  =  3«( LAT*I2+J)-2 

00001660 

0105 

JJ3  =  3*(LAT*I2+J)-5 

00001870 

0106 

JJ4  =  3*{LAT*I+J)-2 

00001880 

0107 

JJ5  =  3*(LAT*I+J)+1 

00001890 

0108 

JJ6  =  3>!'(LAT*I1  +  J)  +  1 

00001900 

0109 

JJ7  =  3#(LAT#12+J)+1 

00001910 

0110 

JJ8  =  34(LAT«Il+J)-5 

00001920 

0111 

JJ9  =  3*{LAT=S'I  +  J)-5 

00001930 

0112 

JJIO  =3*(LAT#Il+J)-8 

00001940 

0113 

JJll  =3>!‘tLAT<‘<  I  +  l)+J)-2 

00001950 

0114 

JJ12  =3*(LAT*Il+J)+4 

00001960 

0115 

JJ13  =3’<'(LAT«(  I-3)+J>-2 

00001970 

0116 

C 

JQl  =  JJl+1 

00001980 

00001990 

0117 

JQ2  =  JJl+2 

00002000 

0118 

C 

DO  102  M=l,  NLAY 

00002010 

00002020 

0119 

IF(M.EQ.1.AND,J.GT.INF( 1) )  GO  TO  102 

00002030 

0120 

IF(M,EQ.l)  GO  TO  192 

00002040 

0121 

IF(J.LE.INF(M-1) .OR.J.GT.INF(M) )  GO  TO  102 

00002050 

0122 

192  IF(J.EQ.l)  GO  TO  200 

00002060 

0123 

IF( I.EQ. IMID.AND.J.EQ.JMID)  GO  TO  203 

00002070 

0124 

IF( I.EQ. IMIO+l.AND. J.EQ. JMID)  GO  TO  203 

00002080 

0125 

IF(J.EQ.LAT)  GO  TO  202 

00002090 

0126 

IF( J,EQ.INF(M) )  GO  TO  201 

00002100 

C 

C 

SHOULD  J  EQUAL  NONE  OF  THE  ABOVE,  CONTINUE  ON  BELOW  TO  STATEMENT  193 

00002110 

00002120 

0127 

C 

193  IFI I.EO.l)  GO  TO  194 

00002130 

00002140 

0128 

IF( I.EQ. FSWl. OR, I.EQ. FSW2)  GO  TO  195 

00002150 

0129 

IF(I.LT.FSW2.AND.I.GT.FSW1)  GO  TO  197 

00002160 

0130 

IF( I.EQ. LAW)  GO  TO  198 

00002170 

C 

C 

EQUILIBRIUM  MATRIX  TERMS  FOR  A  SQUARE  MESH,  H1=H2=H3=H 

00002180 

00002190 

c 

00002200 

0131 

A(KJ1,JJ1>  --8.*(C66(M)+C551M) ) 

00002210 

0132 

A{KJ1,JJ2)  =  4.*C66{M) 

00002220 

0133 

A(KJ1,JJ4)  =  4.»^C66(M) 

00002230 

0134 

A(KJ1,JJ6)  =  4.*C55(M) 

00002240 

0135 

A(KJ1,JJ8)  =  4.=«‘C55<M) 

00002250 

0136 

A(KJ1,JJ1+1)  =  -8.*(C26{M)+C45<M) ) 

00002260 

0137 

A(KJ1,JJ2+1)  =  4.*C26(M) 

00002270 

0138 

A(KJ1,JJ4+1)  =  4.=«'C26(M) 

00002280 

0139 

A(KJ1,JJ6+1)  =  4.*C45{M) 

00002290 

0140 

A(KJ1,JJ8+1)  =  4,*C45(M) 

00002300 

0141 

c 

C  =  C36{M)+C45(M) 

00002310 

00002320 
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0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 


0178 

0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 


C 

A(KJltJJ3+2)  =  C 
A(KJl,JJ5+2)  =  C 
A(KJlrJJ7+2)  =  -C 
A(KJl,JJ9+2)  =  -C 
C 

X(JJl)  =  0. 

C 

A(KQItJJI)  =  -8,*{C26(M}+C45(M) ) 

A(KQ1,JJ2)  =  4,=^C26(M) 

A(KQ1,JJ4)  =  4.>^C26(M) 

A(KQ1,JJ6)  =  4.*C45(M) 

A(KQ1,JJ8)  =  4.>S‘C45(M) 

AIKQlrJJl  +  l)  =  -8.>S'(C22(M)+C44<  M)  ) 

A(KQ1,JJ2+1)  =  4,’J'C22(M) 

A(KQ1,JJ4+1)  =  4,#C22(M) 

A(KQl,JJ6+l)  =  4.#C44(M) 

A{KQ1,JJ8+1)  =  4.*C44(M) 

C 

D  =  C23(M)+C44(M) 

C 

A(K01,JJ3+2)  =  D 
A(KQl»Jj5+2)  =  D 
AIKQl,JJ7+2)  =  -D 
A{KQl,JJ9+2)  =  -D 
C 

X(JOl)  =  0. 

C 

A(KQ2,JJ3)  =  C 
A(KQ2rJJ5)  =  C 
A(KQ2,JJ7)  =  -C 
A<KQ2,JJ9)  =  -C 
A(K02tJJ3+1)  =  D 
A(KQ2,JJ5+1)  =  0 
A(KQ2,JJ7+l)  =  -D 
A{KQ2,JJ9+1)  =  -D 

A(K02,JJ1+2)  =  “8.^(C44(M)+C33(M} J 
A(KQ2tJJ2  +  2)  =  4.=^=C44(M) 

A(KQ2,JJ4+2)  =  4.=i'C44(M) 

A(KQ2,JJ6+2)  =  4,«C33(M) 

A(KQ2,JJ8+2)  =  4,*C33(M) 

C 

X(JQ2)  =  “4.=«'(C13{M)«C2  +  C23(M)*DV  +  2.*C36{  M )  *C4)  =4'HSQ 
GO  TO  102 
C 

C  FREE  SURFACE  MATRIX  TERMS  FOR  1=1  AND  J  NOT  EQUAL  TO  1»  INF  OR  LAT 

194  A(KJltJJl)  =  “3.*C66(M} 

A(KJ1,JJ4}  =  4,’^C66(M) 

A(KJ1,JJ11)  =  -C66{M) 

A(KJ1»JJ1  +  1)  =  -3.=4'C26(M) 

A(KJ1,JJ4+1)  =  4.*C26(M) 

A(KJ1, JJll+1)  =  -C26(M) 

A(KJl,JJ6+2)  =  C36(M) 

A(KJl,JJ8+2)  =  -C36(M) 

C 

A(KQ1,JJ1)  =  -3,=«'C26(M) 

A(KQ1,JJ4)  =4.’!'C26(M) 


00002330 

00002340 

00002350 

00002360 

00002370 

00002380 

00002390 

00002400 

00002410 

00002420 

00002430 

00002440 

00002450 

00002460 

00002470 

00002480 

00002490 

00002500 

00002510 

00002520 

00002530 

00002540 

00002550 

00002560 

00002570 

00002580 

00002590 

00002600 

00002610 

00002620 

00002630 

00002640 

00002650 

00002660 

00002670 

00002680 

00002690 

00002700 

00002710 

00002720 

00002730 

00002740 

00002750 

00002760 

00002770 

00002780 

00002790 

00002800 

00002810 

00002820 

00002830 

00002840 

00002850 

00002860 

00002870 

00002880 

00002890 

00002900 
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0188 

A(KQ1,JJ11)  =  •“C26(M) 

0189 

A(KQ1,JJ1+1)  =  -3.*C22(M) 

0190 

A(KQ1,JJ4+1)  =  4.*C22(M) 

01  1 

A{KQ1»JJH+1J  =  -C22(M) 

01 

A{KQl,JJ6+2)  =  C23<M) 

0193 

A(KQl,JJ8+2)  =  -C23(M) 

0194 

U 

A(KQ2tJJ6)  =  C45IM) 

0?  95 

A(KQ2,JJ8)  =  -C45(M) 

0196 

A(KQ2,JJ6+1)  =  C44(M) 

0197 

A(KQ2,JJ8+1)  =  -C44(M) 

0198 

A<KQ2,JJl  +  2)  =  <-3.*C44{M) 

0199 

A(KQ2,JJ4+2)  =  4.*C44(M} 

0200 

A(KQ2, JJll+2)  =  -C44{M) 

0201 

c 

CYl  =  C12(M)’i'C3  +  C22{M)*BV  +  C26(M)*BU 

0202 

CY2  =  C12(M)#C2  +  C22(M)*DV  +  2  .*C26  (  M ) ’*'64 

0203 

CXYl  =  C16(M)«C3  +  C26(M)’4'BV  +  C66(M)=4'BU 

0204 

CXY2  =  C16{M)«C2  +  C26(M)=«'DV  +  2.*C66(M)*C4 

0205 

c 

X(JJl)  =  -2.*H>!'(CXYi  +  CXY2^Z) 

0206 

X{JQ1}  =  -2.«H*(CY1  +  CY2*Z) 

0207 

X(JQ2)  0. 

0208 

r 

GO  TO  102 

0209 

U 

195  HI  =  H 

0210 

H2  =  FLOAT(K)#H 

0211 

r 

H3  =  H 

0212 

\j 

IF(  I.NE.FSW2)  GO  TO  196 

0213 

HI  =  FLOAT{K)*H 

0214 

c 

H2  -  H 

0215 

196  CONTINUE 

0216 

HH  =  H2/H1 

0217 

HR  =  HH/(1.+HH) 

0218 

HHl  =  H1/H3 

0219 

HH2  =  H2/H3 

0220 

HH3  =  H1«H2 

0221 

HMU  =  HH1<=HH2 

0222 

GO  TO  199 

0223 

197  HI  =  FLOAT{K)=«'H 

0224 

H2  =  HI 

0225 

H3  =  H 

0226 

GO  TO  196 

c 

FREE  SURFACE  MATRIX  TERMS  FOR  I=LAW  AND  J  NOT  EQUAL 

0227 

c 

198  A(KJ1,JJ1)  =  3.*C66(M) 

0228 

A(KJ1,JJ2)  =  -4.’S^C66(M) 

0229 

A(KJ1,JJ13)  =  C66{M) 

0230 

AIKJ1,JJ1+1)  =  3.*C26(M) 

0231 

A(KJ1,JJ2+1)  =  -'4.*C26{M) 

0232 

A(KJ1,JJ13+1)  =  C26(M) 

0233 

A(KJl,JJ6+2)  =  C36(M) 

0234 

A(KJl,JJ8+2)  =  -C36(M) 

C 


00002910 

00002920 

00002930 

00002940 

00002950 

00002960 

00002970 

00002980 

00002990 

00003000 

00003010 

00003020 

00003030 

00003040 

00003050 

00003060 

00003070 

00003080 

00003090 

00003100 

00003110 

00003120 

00003130 

00003140 

00003150 

00003160 

00003170 

00003180 

00003190 

00003200 

00003210 

00003220 

00003230 

00003240 

00003250 

00003260 

00003270 

00003280 

00003290 

00003300 

00003310 

00003320 

00003330 

00003340 

00003350 

00003360 

00003370 

00003380 

00003390 

00003400 

00003^^10 

00003420 

00003430 

00003440 

00003450 

00003460 

00003470 

00003480 
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0235 

A(KQltJJl)  =  3.#C26(M) 

00003490 

0236 

A(KQ1,JJ2)  =  -4,*C26(M) 

00003500 

0237 

A(KQ1,JJ13)  =  C26{M) 

00003510 

0238 

A(KQ1,JJ1+1)  =  3.«C22(M) 

00003520 

0239 

A(KQ1,JJ2+1)  =  ~4.*C22{M) 

00003530 

0240 

A|KQ1,JJ13+1)  =  C22(M) 

00003540 

0241 

A<KQl,JJ6+2)  =  C23(M) 

00003550 

0242 

A(KQl,JJ8+2)  =  ~C23{M) 

00003560 

C 

00003570 

0243 

A(KQ2,JJ6)  =  C45(M) 

00003580 

0244 

A(KQ2,JJ8)  =  -C45(M) 

00003590 

0245 

A(KQ2,JJ6+1)  =  C44(M) 

00003600 

0246 

A(KQ2,JJ8+1)  =  -C44(M) 

00003610 

0247 

A{KQ2tJJl+2)  =  3.*C44(M) 

00003620 

0248 

A{KQ2,JJ2+2)  =  -4.=«‘C44{M) 

00003630 

0249 

A(KQ2tJJI3+2)  =  C44(M) 

00003640 

C 

00003650 

0250 

CYl  =  C12(M)#C3  +  C22(M)>>=BV  +  C26{M)*BU 

00003660 

0251 

CY2  =  C12(M)#C2  +  C22(M)>S‘DV  +  2.>S'C26(  M)*C4 

00003670 

0252 

CXYl  -  C16(M)=^C3  +  C26(M)>«'BV  +  C66(M)’«‘BU 

00003680 

0253 

CXY2  =  C16(M)#C2  +  C26(M)#0V  +  2.«C66{M)*C4 

00003690 

C 

00003700 

0254 

X(JJl)  =  “2.*H*(CXY1  +  CXY2»2) 

00003710 

0255 

XIJQI)  =  «2.*H*ICY1  +  CY2=<‘Z) 

00003720 

0256 

X{JQ2)  =  0. 

00003730 

0257 

GO  TO  102 

00003740 

C 

00003750 

C  EQUILIBRIUM  MATRIX  TERMS  FDR  A  VARIABLE  MESH,  HI, 

H2  f  H3  INDEPENDENT  00003760 

C 

00003770 

0258 

199  A(KJ1,JJ1J  ==  -2.«(C661M)+HMU>5=C55(M)  ) 

00003780 

0259 

A(KJ1,JJ2)  =  2.*HR’8'C66{M) 

00003790 

0260 

A(KJl,JJ4)  =  2.*C66{M)/( l.+HH) 

00003800 

0261 

A{KJ1,JJ6)  =  HMU*C55(M) 

00003810 

0262 

A(KJ1,JJ8)  =  HMU*C55{M) 

00003820 

0263 

AIKJl, JJl+1)  =^2.*(C26(M)+HMU^C45(M) ) 

00003830 

0264 

A{KJI,JJ2+1)  =  2.#HR«C26(M) 

00003840 

0265 

A(KJ1,JJ4+1)  =  2.*C26(M) /( 1,+HH) 

00003850 

0266 

A{KJ1,JJ6+1)  =  HMU=<'C45{M) 

00003860 

0267 

A(KJ1,JJ8+1)  =  HMU=«'C45(M) 

00003870 

C 

00003880 

0268 

C  =  HH1*HR*(C36(M)+C45(M) )/2. 

00003890 

C 

00003900 

0269 

A(KJl,JJ3+2)  =  C 

00003910 

0270 

A(KJl,JJ5+2)  =  C 

00003920 

0271 

A(K-J1,  JJ7+2)  =  -C 

00003930 

0272 

A(KJl,JJ9+2)  =  -C 

00003940 

C 

00003950 

0273 

A(KQ1,JJ1)  =  ~2.*(C26(M)+HMU*C45(M) ) 

00003960 

0274 

A<KQ1,JJ2)  =  2.=«'HR^C26(M) 

00003970 

0275 

AIKQ1,JJ4)  =  2.*C26{M)/ ( l.+HH) 

00003980 

0276 

A(KQ1,JJ6)  =  HMU’#=C45|M) 

00003990 

0277 

A(KQ1,JJ8)  =  HMU*C45(M) 

00004000 

0278 

A(KQ1,JJ1  +  1)  =  -2.=«' (C22tM)+HMU’S'C44<M)  ) 

00004010 

0279 

A(KQ1,JJ2+1)  =  2.=<'HR*C22(M) 

00004020 

0280 

A(KQ1,JJ4+1)  =  2.>>'C22(M)/{  l.+HH) 

00004030 

0281 

A{KQl,JJ6+l)  =  HMU*C44(M) 

00004040 

0282 

A(KQ1,JJ8+1)  =  HMU*C44(M) 

00004050 

C 

00004060 
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0283 

c 

D  =  HH1=«'HR*IC23{M)+C44{M)  )/2. 

0C304070 

00004080 

0284 

A(KQl,JJ3+2)  =  D 

00004090 

0285 

A(KQltJJ5+2>  =  D 

00004100 

0286 

A(KQl,JJ7+2)  =  «D 

00004110 

0287 

c 

A{KQl,JJ9+2)  =  -D 

00004120 

00004130 

0288 

A(KQ2,JJ3)  =  C 

00004140 

0289 

A(KQ2,JJ5)  =  C 

00004150 

0290 

A(KQ2tJJ7)  =  -C 

00004160 

0291 

A(KQ2,JJ9)  =  -C 

00004170 

0292 

A(KQ2,JJ3+1)  =  D 

00004180 

0293 

A(KQ2,JJ5+1)  =  D 

00004190 

0294 

A(KQ2,JJ7+1)  =  -D 

00004200 

0295 

A<KQ2tJJ9+l}  '=  “D 

00004210 

0296 

A{KQ2,JJl+2)  =  -2,*(C44(M)+HMU*C33(M) ) 

00004220 

0297 

A(KQ2,JJ2+2)  =  2.*HR*C44(M) 

00004230 

0298 

A(KQ2»JJ4+2)  =  2.*C44(M)/( l.+HH) 

00004240 

0299 

A(KQ2,JJ6+2)  =  HMU*C33(M) 

00004250 

0300 

c 

A{KQ2,JJ8+2)  =  HMU*C33(M) 

00004260 

00004270 

0301 

X(JJl)  =  0. 

00004280 

0302 

X(JQ1>  -  0. 

00004290 

0303 

X(JQ2)  -  -HH3<'(C13(  M)*C2  +  C23(M)*DV  +  2  .  *C36  ( M )  *C4 ) 

00004300 

0304 

c 

GO  TO  102 

00004310 

00004320 

0305 

200  IF(I.EO.l)  GO  TO  210 

00004330 

0306 

c 

IF{ I.EQ.LAW)  GO  TO  211 

00004340 

00004350 

c 

c 

FREE  SURFACE  MATRIX  TERMS  FOR  I  BETWEEN  1  AND  LAW  AND  J=l. 

00004360 

00004370 

0307 

A(KJ1,JJ1)  =  -3.*C55IM) 

00004380 

0308 

A{KJ1,JJ6)  =  4.*C55iM) 

00004390 

0309 

c 

A(KJ1,JJ121  =  -C55(M) 

00004400 

00004410 

0310 

A(KJ1,JJ1+1)  =  -3.*C45IM) 

00004420 

0311 

A(KJ1,JJ6+1)  =  4.*C45{M) 

00004430 

0312 

c 

A(KJ1,JJ12+1)  =  -C45(M) 

00004440 

00004450 

0313 

A(KQ1,JJ1)  =  ~3.*C45(M) 

00004460 

0314 

A(KQ1,JJ6)  =  4.*C45(M) 

00004470 

0315 

c 

A(KQ1,JJ12)  =  -C45(M) 

00004480 

00004490 

0316 

A(KQ1,JJ1  +  1)  =  -3.’^‘C44IM) 

00004500 

0317 

A(KQ1,JJ6+1)  =  4.#C44(M) 

00004510 

0318 

c 

A(KQ1, JJ12+1)  =  -C44(M) 

00004520 

00004530 

0319 

A(KQ2,JJl+2)  =  -3,*C33(M) 

00004540 

0320 

A(K02tJJ6+2)  =  4.*C33{M) 

00004550 

0321 

c 

A(KQ2, JJ12+2)  =  -C33IM) 

00004560 

00004570 

0322 

CZl  =  ClSIMI^kCS  +  C23<M)«BV  +  C36(M)*BU 

00004580 

0323 

c 

CZ2  =  C13(M)*C2  +  C23(M)*DV  +  2.*C36(M)*C4 

00004590 

00004600 

0324 

X(JJl)  =  0. 

00004610 

0325 

X(JQl)  =  0. 

00004620 

0326 

c 

XIJQ2)  =  -2.«H#(CZ1  +  CZ2*Z) 

00004630 

00004640 
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0327 

IF{ I.EQ.FSWl)  GO  TO  206 

00004650 

0328 

IF( I.EQ.FSW2)  GO  TO  206 

00004660 

0329 

IF( I.GT.FSWl.AND, I.LT,FSW2)  GO  TO 

209 

00004670 

C 

00004680 

C 

IF 

I  IS  BETWEEN  1  AND  FSWl  OR  BETWEEN 

FSW2  AND  LAW,  CONTINUE  BELOW 

00004690 

C 

00004700 

0330 

A(KJl,JJ2+2)  =  -‘C45{M) 

00004710 

0331 

A{KJl,JJ4+2}  =  C45(M} 

00004720 

C 

00004730 

0332 

A(K01,JJ2+2)  =  -C44(M,) 

00004740 

0333 

A{KQl,JJ4+2)  =  C44(M} 

00004750 

C 

00004760 

0334 

A(K02rJJ2)  =  -C36(M) 

00004770 

0335 

A(K02,JJ4)  =  C36(M) 

00004780 

0336 

A(KQ2,JJ2+1)  =  ~C23(M) 

00004790 

0337 

A(KQ2»JJ4+1)  =  C23tM) 

00004800 

0338 

GO  TO  102 

00004810 

C 

00004820 

C 

CASE  WHERE  I=FSW1  OR  FSW2  AND  J=1 

00004830 

C 

00004840 

0339 

206 

XK  =  FLOATIK) 

00004850 

0340 

D1  ^  2.#(XK-1,)/XK 

00004860 

0341 

D2  =  2.«XK/( XK+1. ) 

00004870 

0342 

03  =  2./( (XK+1.)#XK) 

00004880 

C 

00004890 

0343 

IF( I.E0.FSW2)  GO  TO  207 

00004900 

C 

00004910 

0344 

A(KJl,JJl+2)  =  D1^C45(M) 

00004920 

0345 

A(KJl»JJ2+2)  =  -D2=!'C45(M) 

00004930 

0346 

A(KJltJJ4+2)  =  D3*C45{M) 

00004940 

C 

00004950 

0347 

A(KQl,JJl+2)  =  D1*C44(M) 

00004960 

0348 

A{KQl,JJ2+2)  =  “D2*C44(M) 

00004970 

0349 

A(KQl,JJ4+2)  =  D3«C441M) 

00004980 

C 

00004990 

0350 

A(KQ2,JJ1)  =  D1>!'C36(M} 

00005000 

0351 

A(KQ2,JJ2)  =  -D24‘C36(M) 

00005010 

0352 

A(KQ2,JJ4)  =  D3’S=C36(M) 

00005020 

C 

00005030 

0353 

A(KQ2tJJl  +  l)  =  D1=«‘C23(M) 

00005040 

0354 

A(KQ2,JJ2+1)  =  -D2=^=C23(M) 

00005050 

0355 

A{KQ2,JJ4+1)  =  D3«C23(M) 

00005060 

0356 

GO  TO  102 

00005070 

C 

00005080 

0357 

207 

A(KJl,JJl  +  2)  =  -D1=«^C45{M) 

00005090 

0358 

A(KJl,JJ2  +  2)  =  -D3’i=C45(M) 

00005100 

0359 

AIKJl,Jj4+2)  =  D2^C45(M) 

00005110 

C 

00005120 

0360 

A(KQl,JJl+2)  =  -D1«C44(M) 

00005130 

0361 

A{KQl,JJ2  +  2)  =  -D3>}'C44(M) 

00005140 

0362 

A{KQl,Jj4+2)  =  D2*C44(M) 

00005150 

c 

00005160 

0363 

A{KQ2,JJ1)  =  -D1^C36(M) 

00005170 

0364 

A(KQ2,JJ2)  =  -D3*C36(M) 

00005180 

0365 

A(KQ2,JJ4)  =  D2’}'C36{M) 

00005190 

c 

00005200 

0366 

A(KQ2,JJ1+1)  =  “D1=4'C23(M) 

00005210 

0367 

A(KQ2,JJ2+1)  =  -D3*C23(M) 

00005220 
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0368 

A(KQ2,JJ4+1)  =  D2>«'C23(M) 

00005230 

0369 

GO  TO  102 

00005240 

C 

00005250 

C 

CASE 

:  WHERE  I  IS  BETWEEN  FSWl  AND  FSW2  AND  J=1 

00005260 

C 

00005270 

0370 

209 

XK  =  FLOATIK) 

00005280 

0371 

A(KJl,JJ2+2)  =  -C45(M)/XK 

00005290 

0372 

A(KJl,JJ4+2)  =  C45(M)/XK 

00005300 

C 

00005310 

0373 

A(KQl,JJ2+2)  =  -C44(M)/XK 

00005320 

0374 

A(KQUJJ4+2)  =  C44(M)/XK 

00005330 

C 

00005340 

0375 

A(KQ2tJJ2)  =  -C36(M)/XK 

00005350 

0376 

A(KQ2,JJ4I  =  C36(M)/XK 

00005360 

C 

00005370 

0377 

A(KQ2,JJ2+1)  =  *-C23(M)/XK 

00005380 

0378 

A(KQ2tJJ4+l)  =  C23(M)/XK 

00005390 

0379 

60  TO  102 

00005400 

C 

00005410 

C 

FREE  SURFACE  MATRIX  TERMS  FOR  I=J=1 

00005420 

C 

00005430 

0380 

210 

A(KJ1,JJ1)  =  -3.*C66(M) 

00005440 

0381 

A(KJ1,JJ4)  =  4.=«=C66(M) 

00005450 

0382 

AIKJltJJll)  =  -C66IM) 

00005460 

0383 

AIKJltJJl+1)  =  ->3.*C26(M) 

00005470 

0384 

A{KJ1,JJ4+1J  =  4.*C26(M) 

00005480 

0385 

A(KJ1tJJ11+1)  =  «C26(M) 

00005490 

0386 

A(KJl,JJl+2)  =  -3.*C36(M) 

00005500 

0387 

A{KJl,JJ6+2)  4.^C36{M) 

00005510 

0388 

AIKJlt JJ12+2)  =  -C36(M) 

00005520 

C 

00005530 

0389 

AtKQltJJl)  =  -3,#C26(M) 

00005540 

0390 

A(KQltJJ4)  =  4.’}cC26(M) 

00005550 

0391 

A(KQ1,JJ11)  =  -C26(M) 

00005560 

0392 

A{KQ1,JJ1  +  1J  =  “3.=^C22(M) 

00005570 

0393 

A(KQ1,JJ4+1)  =  4.«C22(M) 

00005580 

0394 

A(KQ1,JJ11+1)  =  “C22(M) 

00005590 

0395 

A(KQltJJl+2)  =  “3,*C23(M) 

00005600 

0396 

A(KQl,JJ6+2)  =  4.*C23(M) 

00005610 

0397 

A(KQ1,  JJ12+2)  =  >-C23(M) 

00005620 

C 

00005630 

0398 

A{KQ2,JJ1)  =  -‘3.*C45(M) 

00005640 

0399 

A(K02tJJ6)  =  4,*C45(M) 

00005650 

0400 

A(KQ2,JJ12)  =  -C45(M) 

00005660 

0401 

A(KQ2,JJ1+1)  =  -3.#C44(M) 

00005670 

0402 

A<KQ2tJJ6+l)  =  4,*C44(M) 

00005680 

0403 

A(KQ2,JJ12+1}  =  -C44(M) 

00005690 

0404 

A(KQ2tJJl+2)  =  -3.*C44{M) 

00005700 

0405 

A(KQ2»JJ4+2)  =  4,*C44(M) 

00005710 

0406 

A{KQ2, JJll+2)  =  ~C44(M) 

00005720 

C 

00005730 

0407 

CYl  =  C12(M)*C3  +  C22(M)*BV  +  C26(M)=«‘BU 

00005740 

0408 

CY2  =  C12(M)«C2  +  C22{M)*DV  +  2.=f'C26{  M)*C4 

00005750 

0409 

CXYl  =  C16(M)*C3  +  C26(M)#BV  +  C66(M)*BU 

00005760 

0410 

CXY2  =  C16(M)*C2  +  C26(M)#DV  +  2,*C66(M)*C4 

00005770 

C 

00005780 

0411 

X(JJl)  =  -2,«H«(CXY1  +  CXY2*Z) 

00005790 

0412 

X(JQl)  ==  -2.>tcH*(CYl  +  CY2*Z) 

00005800 
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0413 

X(JQ2)  =  0. 

00005810 

0414 

GO  TO  102 

00005820 

C 

00005830 

C 

FREE  SURFACE  MATRIX  TERMS  FOR  J=1  AND  I=LAW 

00005-840 

C 

00005850 

0415 

211 

A(KJ1»JJ1)  -  3.#C66(M) 

00005860 

0416 

A(KJ1,JJ2)  =  -4.#C66(M) 

00005870 

0417 

A(KJ1,JJ13)  =  C66{M) 

00005880 

0418 

A(KJ1,JJ1+1)  =  3.’!'C26IM) 

00005890 

0419 

A(KJ1,JJ2+1)  =  -4.’;‘C26(M) 

00005900 

0420 

A(KJ1, JJ13+1)  =  C26(M) 

00005910 

0421 

A<KJl,JJl+2)  =  ~3.*C36(M) 

00005920 

0422 

A(KJl,JJ6+2)  =  4.*C36(M) 

00005930 

0423 

A(KJl,JJ12+2)  =  -C36(M) 

00005940 

C 

00005950 

0424 

A(KQlfJJl)  =  3.*C26(M) 

00005960 

0425 

A(KQ1,JJ2)  =  -“4,*C26(M) 

00005970 

0426 

A(KQltJJ13)  =  C26(M) 

00005980 

0427 

A{KQ1,JJ1+1)  =  3,*C22(M) 

00005990 

0428 

A(KQ1,JJ2  +  1)  =  -4.’«'C22(M) 

00006000 

0429 

A(K01, JJ13+1)  =  C22(M) 

00006010 

0430 

A(KQl,JJl  +  2)  =  -3.’«'C23(M) 

00006020 

0431 

A(KQl,JJ6+2)  =  4.«C23(M) 

00006030 

0432 

A(K01»JJ12+2)  =  -C23(M) 

00006040 

C 

00006050 

0433 

A(KQ2»JJ1)  -  -3.4C45IM) 

00006060 

0434 

A(KQ2,JJ6)  =  4,*C45IM) 

00006070 

0435 

A(KQ2,JJ12)  =  -C45(M) 

00006080 

0436 

A(KQ2,JJl+l)  =  -3.*C44(M) 

00006090 

0437 

A(KQ2rJJ6+l)  =  4.*C44(M) 

00006100 

0438 

At KQ2, JJ12+1 )  =  -C44(M) 

00006110 

0439 

A(KQ2tJJ1+2)  =  3.*C44{M) 

00006120 

0440 

A(KQ2,JJ2+2)  =  -4,*C44(M) 

00006130 

0441 

A(KQ2,JJ13+2)  =  C44(M) 

00006140 

C 

00006150 

0442 

CYl  =  C12(M)*C3  +  C22(M)#BV  +  C26(M)*BU  . 

00006160 

0443 

CY2  =  C12(M)*C2  +  C22<M)*DV  +  2.  =*^026  (  M )  =4^04 

00006170 

0444 

CXYl  =  C16(M)*C3  +  C26(M)=«=BV  +  C66{M)*BU 

00006180 

0445 

CXY2  =  C16(M)»!‘C2  +  C26(M)*DV  +  2.*C66(M)*C4 

00006190 

C 

00006200 

0446 

XtJJl)  =  -2.«H*(CXY1  +  CXY2=«'Z) 

00006210 

0447 

X(JQl)  =  -2.*H*{CY1  +  CY2*Z) 

00006220 

0448 

X{J02)  =  0. 

00006230 

0449 

GO  TO  102 

00006240 

C 

00006250 

0450 

201 

P  =  M+1 

00006260 

0451 

IF( I.EQ. 1)  GO  TO  220 

00006270 

0452 

IF( I.EQ. FSWl)  GO  TO  221 

00006280 

0453 

IF( I,LT.FSW2.AND.I.GT.FSWl)  GO  TO  222 

00006290 

0454 

IF{  I.EQ.FSW2)  GO  TO  221 

00006300 

0455 

IF( I.EQ. LAW)  GO  TO  223 

00006310 

C 

00006320 

C 

MATRIX  TERMS  AT  INTERFACE  FOR  I  BETWEEN  1  AND  FSWl  OR  FSW2  AND  LAW  00006330 

C 

00006340 

0456 

A(KJ1,JJ1)  =  3.*{C55(M)+C55(P) ) 

00006350 

0457 

A|KJ1,JJ6)  =  ->4.>«'C55(  P) 

00006360 

0458 

A(KJ1,J.J8)  =  -4,’i'C55(M) 

00006370 

0459 

AtKJl,JJ10)  =  C55(M) 

00006380 
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0460 

A(KJ1,JJ12)  =  C55(P) 

00006390 

C 

00006400 

0461 

A(KJ1,JJ1+1)  =  3.*(C45(M)+C45IP) ) 

00006410 

0462 

A(KJ1,JJ6+1)  =  -4.’«'C45(P) 

00006420 

0463 

A(KJltJJ8+l)  =  -4.*C45<M) 

00006430 

0464 

A(KJ1, JJlO+l)  =  C45(M) 

00006440 

0465 

A(KJ1,JJ12+1)  =  C45{P) 

00006450 

C 

00006460 

0466 

A(KJl,JJ2+2)  =  C45(P)-C45(M) 

00006470 

0467 

A(KJl,JJ4+2)  =  C45(M)-C45( P) 

00006480 

C 

00006490 

0468 

A(KQ1,JJ1)  =  3.>«‘(C45(M)+C45(P)  ) 

00006500 

0469 

A(KQ1,JJ6)  =  -4.#C45(P) 

00006510 

0470 

A(KQ1,JJ8)  =  ~4.*C45(M) 

00006520 

0471 

A(KQ1,JJ10)  =  C45(M) 

00006630 

0472 

A(KQ1,JJ12)  =  C45(P) 

00006540 

C 

00006550 

0473 

AIKQ1,JJ1+1)  =  3.’i'(C44(MI  +  C44(  P)  ) 

00006560 

0474 

A(KQ1,JJ6+1)  =  -4.«C44(P) 

00006570 

0475 

A(KQ1,JJ8+1)  =  -4.*C44<M) 

00006580 

0476 

AIKQ1,JJ10+1)  =  C44IM) 

00006590 

0477 

A(KQ1, JJ12+1)  =  C44(P) 

00006600 

C 

00006610 

0478 

A(K01,JJ2+2)  =  C44( P)-C44(M) 

00006620 

0479 

A{KQl,JJ4+2)  =  C44(M)~C44(P) 

00006630 

c 

00006640 

0480 

A{KQ2rJJ2)  =  C36 ( P ) -C36 ( M ) 

00006650 

0481 

A(KQ2,JJ4)  -  C36(M)-C36(P) 

00006660 

c- 

00006670 

0482 

A{KQ2,JJ2+1)  =  C23(  P)-’C23(M) 

00006680 

0483 

AIKQ2,JJ4+1)  =  C23(M)-C23<P> 

00006690 

c 

00006700 

0484 

AlKQ2,JJl+2)  =  3.*{C33(M)+C33(P) ) 

00006710 

0485 

A(KQ2,JJ6+2)  =  -4.*C33(P) 

00006720 

0486 

A(KQ2tJJ8+2)  =  -4.*C33(M) 

00006730 

0487 

A(KQ2, JJlO+2)  =  C33(M) 

00006740 

0488 

A(KQ2,JJ12+2)  =  C33IP) 

00006750 

c 

00006760 

0489 

CZl  =  (C13I  P)-C13(M)  »*C3  +  (C23(  P)~C23(M)  )*BV  +  (  C36  (  P  ) -C36  ( M  )  )  >i'BU00006770 

0490 

CZ2  =  (C13{P)-C13(M) )*C2+(C23{ P)-C23(M) )*DV+2.*( C36{ P )-C36( M J )*C4 

00006780 

c 

00006790 

0491 

XIJJl)  =  0. 

00006800 

0492 

X(JQl)  =  0. 

00006810 

0493 

X(J02)  =  2.*H*{CZ1  +  CZ2*Z) 

00006820 

0494 

GO  TO  102 

00006830 

c 

00006840 

c 

FREE  SURFACE  MATRIX  TERMS  AT  ANY  INTERFACE  WHERE  1=1  AND  J=INF  OR  AT 

00006850 

c 

THE 

FREE  SURFACE  POINT  1=1,  J=LAT 

00006860 

c 

00006870 

0495 

220 

A(KJltJJl)  =  “3.*C66{M) 

00006880 

0496 

A(KJ1,JJ4)  =  4,*C66(M) 

00006890 

0497 

A(KJ1,JJ11)  =  --C66(M) 

00006900 

c 

00006910 

0498 

A(KJ1,JJ1+1)  =  ~3.#C26(M) 

00006920 

0499 

A(KJ1,JJ4+1)  =  4,#C26(M) 

00006930 

0500 

AUJ1,JJ11+1)  =  -C26(M) 

00006940 

c 

00006950 

0501 

A{KJl,JJl+2)  =  3.*C36(M) 

00006960 

59 


FORTRAN  IV  G1  RELEASE  2.0 


MAIN 


DATE  =  75007 


08/16/07 


0502 

0503 

0504 

0505 

0506 

0507 

0508 

0509 

0510 

0511 

0512 

0513 

0514 

0515 

0516 

0517 

0518 

0519 

0520 

0521 

0522 

0523 

0524 

0525 

0526 

0527 

0528 

0529 


0530 

0531 

0532 

0533 

0534 

0535 

0536 

0537 

0538 

0539 

0540 

0541 

0542 

0543 

0544 

0545 


A(KJl,JJ8+2)  =  ~4.«C36(M) 

A(KJl,JJ10+2)  =  C36(M) 

C 

A(KQ1,JJ1)  -  -3.^C26(M) 

A(K01,JJ4)  =4.’«'C26(M) 

A(KQ1,JJ11)  =  -C26(M) 

C 

A(KQ1,JJ1+1)  =  “3,^C22(M) 

A(KQ1,JJ4+1)  =  4,*C22(M) 

A{KQ1,JJ11+1)  =  -C22(M) 

C 

A(KQl,JJl  +  2)  =  3.:i^C23(M) 

A(KQl,JJ8+2)  =  ~4.’^C23(M) 

A(KQl,JJ10+2)  =  C23(M) 

C 

A(KQ2,JJ1)  =  3.=«'C45(M) 

A(KQ2,JJ8}  =  -4.»!=C45<M) 

A(KQ2tJJ10)  =  C45(M) 

C 

A(KQ2,JJ1+1)  =  3,=i‘C44(M) 

A(KQ2,JJ8+1)  =  -4.=(‘C44(M) 

A(KQ2,JJ10+1)  =  C44(M) 

C 

A(KQ2,JJl+2)  =  -3.*C44(M) 

A(KQ2,JJ4+2)  =  4,«C44(M) 

A(KQ2,JJll+2)  =  -C44(M) 

C 

CYl  =  C12(M)=i‘C3  +  C22(M)=«‘BV  +  C26(M)=«'6U 
CY2  =  C12{M)^C2  +  e22(M)*DV  +  2  •=5‘C2  6  (  M ) ’S'CA 
CXYl  =  C16(M)’!'C3  +  C26(M)*BV  +  C66{M)*BU 
CXY2  =  C16(M)=«=C2  +  C26{M)’!'DV  +  2. ’J'C  66  (  M )  >l‘C4 
C 

X(JJl)  =  -2.*HX=(CXY1  +  CXY2*Z) 

X(JQl)  =  -2.^H=!'{CY1  +  CY2*Z) 

XIJQ2)  =  0. 

GO  TO  102 
C 

C  MATRIX  TERMS  AT  THE  INTERFACE  FOR  J=INF  AND  I=FSW1  OR  I=FSW2 
C 

221  XK  =  FLOAT(K) 

Dl  =  (XK-1.)/XK 
D2  =  XK/(XK+1,) 

D3  =  l./( (XK+1. )^XK) 

C 

A(KJ1»JJ1)  =  3.’!‘{C55(M)+C55{  P)  ) 

A(KJ1,JJ6)  =  -4.=<'C55(P) 

A(KJ1,JJ8)  =  -4.*C55(M) 

A(KJ1,JJ10)  =  C55(M) 

A(KJ1, JJ12)  =  C55( P) 

C 

A(KJUJJ1+1)  =  3.=S'(C45(M)+C45(  P)  ) 

A(KJ1,JJ6+1)  =  -4,’!'C45(P) 

A(KJ1,JJ8+1)  =  -4.«C45(M) 

A(KJ1,JJ10+1)  =  C45(M) 

A(KJ1, JJ12+1)  =  C45(P) 

C 

A(KQItJJI)  =  3.*(C45(M)  +  C45(P)) 

A{KQ1,JJ6)  =  -4.''!'C45  (  P) 


00006970 

00006980 

00006990 

00007000 

00007010 

00007020 

00007030 

00007040 

00007050 

00007060 

00007070 

00007080 

00007090 

00007100 

00007110 

00007120 

00007130 

00007140 

00007150 

00007160 

00007170 

00007180 

00007190 

00007200 

00007210 

00007220 

00007230 

00007240 

00007250 

00007260 

00007270 

00007280 

00007290 

00007300 

00007310 

00007320 

00007330 

00007340 

00007350 

00007360 

00007370 

00007380 

00007390 

00007400 

00007410 

00007420 

00007430 

00007440 

00007450 

00007460 

00007470 

00007480 

00007490 

00007500 

00007510 

00007520 

00007530 

00007540 
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0546 

AIKQUJJS)  =  -4.’i'C45(M) 

00007550 

0547 

A(KQ1,JJ10)  =  C45(M) 

00007560 

0548 

A(KQ1,JJ12)  =  C45(P) 

00007570 

C 

00007580 

0549 

A|KQ1,JJ1+1)  =  3.*{C44{M)+C44I 

[P)  ) 

00007590 

0550 

A(KQ1,JJ6+1)  =  -4,>:'C44(P) 

00007600 

0551 

A{KQ1,JJ8+1)  =  -4.*C44(M) 

00007610 

0552 

A(KQ1, JJlO+1)  =  C44(M) 

00007620 

0553 

A(KQ1,JJ12+1)  =  C44(P) 

00007630 

C 

00007640 

0554 

A(KQ2,JJl  +  2)  ==  3.*(C33{M)+C33( 

[P)  ) 

00007650 

0555 

A{KQ2,JJ6+2)  =  -4.=<«C33|P) 

00007660 

0556 

A(KQ2,JJ8+2)  =  -4.*C33{M) 

00007670 

0557 

A{KQ2, JJlO+2)  =  C33{M) 

00007680 

0558 

A(KQ2,JJ12+2)  =  C33(P) 

00007690 

C 

00007700 

0559 

CZl  =  (C13(P)-C13(M) )*C3  +  (C23(P) 

-C23(M))*BV  +  (C36(P)-C36(M) )*BU00007710 

0560 

CZ2  =  (C13(  P}-*C13(M}  P)-C23(M)  )^DV+2.>!'(C36I  P)-C36(M)  )>!'C4  00007720 

C 

00007730 

0561 

X(JJl)  =  0. 

00007740 

0562 

XlJOl)  =  0. 

00007750 

0563 

X(JQ2)  =  2.«H=!‘{CZ1  +  CZ2*Z) 

00007760 

C 

00007770 

0564 

C=C45(M)^C45(P} 

•  00007780 

0565 

D=C44(M )-C44( P) 

00007790 

0566 

E=C23(M)-C23(P) 

00007800 

0567 

CC=C36(M)-C36(P) 

00007810 

C 

00007820 

0568 

IF( I,EQ.FSW2)  GO  TO  227 

00007830 

C 

00007840 

0569 

A(KJl,JJl  +  2)  =  Z.^J'Dl^C 

00007850 

0570 

A(KJl,JJ2+2)  =  -2.*D2=4=C 

00007860 

0571 

A(KJl,JJ4+2)  =  2.*D3*C 

00007870 

C 

00007880 

0572 

A(KQl,JJl  +  2)  =  Z.^Dl’J'D 

00007890 

0573 

A{KQl,JJ2+2)  =  -2,*D2*D 

00007900 

0574 

A(KQl,JJ4+2)  =  2. =<‘03*0 

00007910 

C 

00007920 

0575 

A(KQ2tJJl)  =  2.*D1*CC 

00007930 

0576 

A{KQ2»JJ2)  =  -2. *02=5=00 

00007940 

0577 

A(KQ2,JJ4)  =  2. *03*00 

00007950 

C 

00007960 

0578 

A(KQ2»JJ1  +  1)  ==  2.*D1=!‘E 

00007970 

0579 

A(KQ2,JJ2+1)  =  -2.*D2*E 

00007980 

0580 

A(KQ2,JJ4+1)  =  2.*D3*E 

00007990 

0581 

GO  TO  102 

00008000 

C 

00008010 

05  82 

227 

A(KJl,JJl+2)  =  “2. *01*0 

00008020 

0583 

A(KJltJJ2+2)  =  -2. *03*0 

00008030 

0584 

A(KJl»JJ4+2)  =  2. *02*0 

00008040 

C 

00008050 

0585 

A(KQl,JJl+2)  =  -2. *01*0 

00008060 

0586 

A(KQl»JJ2+2)  =  -2. *03*0 

00008070 

0587 

A(KQl,JJ4+2)  =  2. *02*0 

00008080 

C 

00008090 

0588 

A(KQ2,JJ1)  =  -2. *01*00 

00008100 

0589 

AUQ2,JJ2)  =  -2. *03*00 

00008110 

0590 

A(KQ2,JJ4)  =  2. *02*00 

00008120 
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0591 

0592 

0593 

0594 


0595 

0596 

0597 

0598 

0599 

0600 

0601 

0602 

0603 

0604 

0605 

0606 

0607 

0608 

0609 

0610 

0611 

0612 

0613 

0614 

0615 

0616 

0617 

0618 

0619 

0620 

0621 

0622 

0623 

0624 

0625 

0626 

0627 

0628 

0629 

0630 

0631 

0632 

0633 

0634 


A{KQ2,JJ1  +  1)  =  -2.^D1»!'E 
A(KQ2,JJ2+1)  =  '-2.*D3*E 
A{KQ2,JJ4+1)  =  2.*D2>!‘E 
GO  TO  102 


MATRIX  TERMS  AT  AN  INTERFACE  FOR  J^^INF  AND  I  BETWEEN  FSWl  AND  FSW2 


222  XK  =  FLOAT{K) 

A(KJI,JJ1)  =  3.*(C55(M)+C55( P) ) 
A{KJ1,JJ6)  =  -4.*C55{P) 

A(KJ1tJJ8)  =  -4,*C55(M) 
A(KJl,JJ10)  =  C55(M) 

A(KJ1,JJ12)  =  C55(P) 

A(KJ1,JJ1+1)  =  3,#(C45(M)+C45IP) ) 
A(KJ1,JJ6+1)  =  -4,*C45{P) 
A{KJltJJ8+l)  =  -4.^C45(M) 

A(KJ1, JJlO+1)  =  C45(M) 
A(KJ1,JJ12+1)  =  C45(P) 

A(KJl,JJ2+2)  =  (C45( P)~C45(M) )/XK 
A(KJl,JJ4+2)  =  {C45{M)-C45|P) )/XK 

A(K01,JJ1)  =  3.4(C45(M)+C45( P) ) 
A(KQ1,JJ6)  =  -4.>CcC45{P) 

A(KQ1»JJ8)  =  -4.4C45(M) 
A(KQ1,JJ10)  =  C45(M) 

A{KQ1,JJ12)  =  C45<P) 

A(KQ1,JJ1+1)  =  3.*(C44(M)+C44( P) ) 
A(KQltJJ6+l)  =  -4.>!‘C44(P) 
A(KQ1,JJ8+1)  =  -4.*C44(M) 
A(KQ1,JJ10+1)  =  C44(M) 

A(KQ1, JJ12+1)  =  C44(P) 

A(KQl,JJ2+2)  =  (C44(P)-C44{M) )/XK 
A(KQl,JJ4+2)  =  (C44(M)-C44(P) )/XK 

A(KQ2,JJ2)  =  (C36(P)-C36(M) )/XK 
A(K02,JJ4)  =  {C36{M)-C36(P) )/XK 


A<KQ2, JJ2+1) 
A(KQ2,JJ4+1) 


(C23( P)-C23(M) ) /XK 
{C23(M)-C23{P) )/XK 


A(KQ2,JJl+2)  =  3.#(C33(M)+C33{P) ) 
A(KQ2,JJ6+2)  =  -4,>S'C33(P) 

A(KQ2,JJ8+2)  =  -4,’!'C33(M) 

A{KQ2, JJlO+2)  =  C33{M) 

A(KQ2,JJ12+2)  =  C33(P) 

X{JQ1)  =  0. 

CZl  =  (C13(P)'-C13(M)  )*C3  +  {C23(P)*C23( 
CZ2  =  (C13(P)-C13(M) )*C2+(C23(P)-C23(M) 

X( JJl )  =  0. 

X{JQ2)  =  2.*H«{CZ1  +  CZ2*Z) 

GO  TO  102 


00008130 

00008140 

00008150 

00008160 

00008170 

00008180 

I  BETWEEN  FSWl  AND  FSW2  00008190 

00008200 

00008210 

00008220 

00008230 

00008240 

00008250 

00008260 

00008270 

00008280 

00008290 

00008300 

00008310 

00008320 

00008330 

00008340 

00008350 

00008360 

00008370 

00008380 

00008390 

00008400 

00008410 

00008420 

00008430 

00008440 

00008450 

00008460 

00008470 

00008480 

00008490 

00008500 

00008510 

00008520 

00008530 

00008540 

00008550 

00008560 

00008570 

00008580 

00008590 

00008600 

00008610 

00008620 

00008630 

00008640 

M))«BV  +  (C36{ P)-C36<M) )^BU00008650 

)^DV+2.*(C36( P)-C36(M) )«C4  00008660 

00008670 

00008680 

00008690 

00008700 
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C 

00008710 

C 

FREE  SURFACE  MATRIX  TERMS  AT  ANY  INTERFACE  WHERE  I=LAW»  J=INF  OR  AT 

00008720 

C 

THE 

FREE  SURFACE  POINT  I=LAW,  J=LAT 

00008730 

c 

00008740 

0635 

223 

AlKJltJJl)  =  3.*C66IM) 

00008750 

0636 

A(KJ1,JJ2)  =  -4.^C66(M) 

00008760 

0637 

A(KJlrJJ13)  =  C66(M) 

00008770 

c 

00008780 

0638 

A(KJ1,JJ1+1)  =  3.«C26(M) 

00008790 

c 

00008860 

0644 

A(KQ1,JJ1)  =  3.*C26{M) 

00008870 

0645 

A(KQ1,JJ2)  =  ~4.«C26(M) 

00008880 

0646 

A(KQ1,JJ13)  =  C26(M) 

00008890 

c 

00008900 

0647 

A<KQ1,JJI+1)  =  3.*C22{M) 

00008910 

0648 

AIKQ1,JJ2+1)  =  ~4.#C22(M) 

00008920 

0649 

A(KQ1,JJ13+1}  =  C22(M) 

00008930 

c 

00008940 

0650 

A(KQUJJl  +  2)  =  3.«C23(M) 

00008950 

0651 

A{K01,JJ8+2)  -  -4.«C23(M) 

00008960 

0652 

A(KQlt JJlO+2)  =  C23(M) 

00008970 

c 

00008980 

0653 

A{KQ2tJJlI  =  3.*C45(M) 

00008990 

0654 

A(KQ2,JJ8)  -  -4.*C45(M) 

00009000 

0655 

A(KQ2,JJ10)  =  C45(M) 

00009010 

c 

00009020 

0656 

A(KQ2,JJ1+1)  =  3.*C44(M) 

00009030 

0657 

A(KQ2tJJ8+1)  =  -4.*C44(M) 

00009040 

0658 

A(KQ2,JJ10+l)  =  C44(M) 

00009050 

c 

00009060 

0659 

A(KQ2tJJl+2)  =  3.«C44(M) 

00009070 

0660 

AtKQ2tJJ2+2)  =  ”4,=i«C44{M) 

00009080 

0661 

A(KQ2,JJ13+2)  =  C44{M) 

00009090 

c 

00009100 

0662 

CYl  =  C12(M)>!'C3  +  C22(M)>!‘BV  +  C26(M)=«'BU 

00009110 

•0663 

CY2  =  C12(M)*C2  +  C22(M)*DV  +  2**C26I  M ) ’(‘C4 

00009120 

0664 

CXYl  =  C16<M)=)«C3  +  C26(M)*BV  +  C66(M)=!‘BU 

00009130 

0665 

CXY2  =  C16(M)^C2  +  C26(M)*DV  +  2.*C66(M)*C4 

00009140 

c 

00009150 

0666 

X(JJl)  =  “2,«H*(CXYi  +  CXY2*Z) 

00009160 

0667 

X(JQl)  =  ~2.*H*(CY1  +  CY2=<^Z) 

00009170 

0668 

X(JQ2)  =  0. 

00009180 

0669 

GO  TO  102 

00009190 

c 

00009200 

c 

MATRIX  TERMS  TO  FIX  THE  RIGID  TRANSLATIONS 

00009210 

c 

00009220 

0670 

203 

A(KJ1,JJ1)  =  1.0 

00009230 

0671 

AIKQltJJl+1)  =  1.0 

00009240 

0672 

A(KQ2tJJ1+2)  =  1.0 

00009250 

c 

00009260 

0673 

X(JJl)  =  0. 

00009270 

0674 

XIJQI)  =  0. 

00009280 
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0675 

0676 

0677 

0678 


0679 

0680 

0681 

0682 

0683 

0684 

0685 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 


0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 


0711 

0712 


X(JQ2)  =  0. 

GO  TO  102 
C 

202  IF( I.EQ.l)  GO  TO  220 
IF( I.EQ.LAW)  GO  TO  223 
C 

C  FREE  SURFACE  MATRIX  TERMS  FOR  I  BETWEEN  1  AND  LAW  AND  J=LAT 
C 

A(KJ1,JJ1)  =  3,*C55(M) 

A(KJltJJ8)  =  -4.’«'C55(M) 

A(KJ1,JJ10)  =  C55(M) 

C 

A(KJ1,JJ1+1)  =  3.*C45(M) 

A(KJ1,JJ8+1)  =  -4.’!‘C45(M) 

A(KJ1, JJlO+1 )  =  C45(M) 

C 

A(K01,JJ1)  =  3.=5'C45{M) 

A(KQ1,JJ8)  =  -4.*C45(M) 

A(KQ1»JJ10)  ==  C45(M) 

C 

A(KQ1,JJ1  +  1)  =  3.>«‘C44(M) 

A(K01,JJ8+1)  =  -4.>i‘C44(M) 

A(KQ1, JJlO+1 )  =  C44(M) 

C 

A(KQ2,JJl  +  2)  =  3.’S'C33(M) 

A(KQ2tJJ8+2)  =  -4,#C331M) 

A(KQ2, JJlO+2)  =  C331M) 

C 

CZl  =  C13IM)*C3  +  C23(M}=i'BV  +  C36<M)’«‘BU 
CZ2  =  C13{M)*C2  +  C23(M)’i‘DV  +  2.«C36(M)*C4 
C 

X(JJl)  =  0. 

X(JOl)  =  0, 

X(JQ2I  =  -Z.^H’KICZI  +  CZ2=^Z) 

C 

IF( I .EQ.FSWl)  GO  TO  231 
IF( I.EQ.FSW2)  GO  TO  231 
IF( I,GT.FSW1.AND.I.LT.FSW2)  GO  TO  234 
C 

C  IF  I  IS  BETWEEN  1  AND  FSWl  OR  BETWEEN  FSW2  AND  LAW,  CONTINUE  BELOW 
C 

A(KJl,JJ2+2)  =  -C45{M) 

A(KJl,JJ4+2)  =  C45(M) 

C 

A(K01,JJ2+2)  =  -C44(M) 

A(KQl,JJ4+2)  =  C44(M) 

C 

A(K02,JJ2)  =  -C36(M) 

A(KQ2,JJ4)  =  C36IM) 

C 

A(KQ2,JJ2+1)  =  -C23(M} 

A(KQ2,JJ4+1)  =  C23(M) 

GO  TO  102 
C 

C  CASE  WHERE  I=FSW1  OR  FSW2  AND  J=LAT 
C 

231  XK  =  FLOAT(K) 

D1  =  2.^MXK^1.)/XK 


00009290 

00009300 

00009310 

00009320 

00009330 

00009340 

00009350 

00009360 

00009370 

00009380 

00009390 

00009400 

00009410 

00009420 

00009430 

00009440 

00009450 

00009460 

00009470 

00009480 

00009490 

00009500 

00009510 

00009520 

00009530 

00009540 

00009550 

00009560 

00009570 

00009580 

00009590 

00009600 

00009610 

00009620 

00009630 

00009640 

00009650 

00009660 

00009670 

00009680 

00009690 

00009700 

00009710 

00009720 

00009730 

00009740 

00009750 

00009760 

00009770 

00009780 

00009790 

00009800 

00009810 

00009820 

00009830 

00009840 

00009850 

00009860 
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0713 

D2  =  2.’S'XK/(XK+1.) 

00009870 

0714 

D3  =  2./( (XK+1.)*XK) 

00009880 

C 

00009890 

0715 

IF( I.EQ.FSW2] 

t  GO  TO  232 

00009900 

C 

00009910 

0716 

AtKJl, JJl+2 J 

=  D1#C45{M) 

00009920 

0717 

A{KJl,JJ2-«-2) 

=  -’D2>i'C45(M) 

00009930 

0718 

A(KJ1, JJ4+2) 

=  D34C45(M) 

00009940 

C 

00009950 

0719 

A(KQ1, JJl+2) 

=  D1*C44(M) 

00009960 

0720 

A(KQl,JJ2+2) 

=  -D2«C44(M) 

00009970 

0721 

AUQl,  JJ4+2) 

=  D3#C44(M) 

00009980 

C 

00009990 

0722 

A(KQ2tJJl)  = 

D1>1‘C36(M) 

00010000 

0723 

A(KQ2,JJ2)  = 

-D2*C36(M) 

00010010 

0724 

AIKQ2,JJ4)  = 

D3#C36(M) 

00010020 

c 

00010030 

0725 

A(KQ2, JJl+1 ) 

=  D1«C23(M) 

00010040 

0726 

A(K02tJJ2+1) 

=  ■'D2*C23(M) 

00010050 

0727 

A(K02, JJ4+1 ) 

-  D3*C23(M) 

00010060 

0728 

GO  TO  102 

00010070 

c 

00010080 

0729 

232 

A(KJ1, JJl+2) 

=  -D1*C45(M) 

00010090 

0730 

A(KJ1, JJ2+2) 

=  ”D3=i'C45(M) 

00010100 

0731 

A(KJl,JJ4+2) 

=  02*C45<M) 

00010110 

c 

00010120 

0732 

AIKQl, JJl+2) 

=  -D1’!'C44(M) 

00010130 

0733 

AIKQl, JJ2+2) 

=  -D3*C44(M) 

00010140 

0734 

A(KQ1 ,JJ4+2) 

=  D2*C44(M) 

00010150 

c 

00010160 

0735 

A(KQ2,JJ1)  = 

-D1>4‘C36(M) 

00010170 

0736 

A(KQ2,JJ2)  = 

-03*C36(M) 

00010180 

0737 

A(KQ2,JJ4)  = 

D2*C36IM) 

00010190 

c 

00010200 

0738 

A(KQ2,JJ1+1) 

=  -D1*C23|M) 

00010210 

0739 

A(KQ2, JJ2+1) 

=  -D3*C23(M) 

00010220 

0740 

A(KQ2,JJ4+1) 

=  D24C23(M) 

00010230 

0741 

GO  TO  102 

00010240 

c 

00010250 

c 

CASE  WHERE  I  IS  BETWEEN  FSWl  AND  FSW2  AND  J=LAT 

00010260 

c 

00010270 

0742 

234 

XK  =  FLOAT(K) 

00010280 

0743 

A(KJl,JJ2+2) 

=  -C45(M)/XK 

00010290 

0744 

A(KJ1, JJ4+2) 

=  C45(M)/XK 

00010300 

c 

00010310 

0745 

A(KQ1, JJ2+2) 

=  -C44(M)/XK 

00010320 

0746 

A(K01 ,JJ4+2) 

=  C44(M)/XK 

00010330 

c 

00010340 

0747 

A(KQ2,JJ2)  = 

-C36(M)/XK 

00010350 

0748 

A(KQ2,JJ4)  = 

C36I M) /XK 

00010360 

c 

00010370 

0749 

A(KQ2, JJ2+1) 

=  -C23(M)/XK 

00010380 

0750 

A(KQ2,JJ4+1) 

=  C231M}/XK 

00010390 

c 

00010400 

0751 

102 

CONTINUE 

00010410 

c 

00010420 

c 

FORM  THE  NDNSYMETRIC  BANDED  MATRIX  AX 

00010430 

c 

00010440 
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0752 

IL  =  KJ1+3*!N0DE-1) 

0753 

C 

IN  =  IL+2 

0754 

DO  103  IK=IL,  IN 

0755 

C 

II  =  IK-IL+1 

0756 

DO  104  JK=1, NBAND 

0757 

JJ  =  IK+JK-IBW-1 

0758 

IF! IK.LE.IBWl)  JJ  =  JK 

0759 

IF! JJ.GT. JQMAX)  GO  TO  105 

0760 

AX! JK,  IK)  =  A! II,JJ) 

0761 

GO  TO  104 

0762 

105 

AX!JK,IK)  =  0,0 

0763 

104 

CONTINUE 

0764 

103 

CONTINUE 

0765 

101 

CONTINUE 

0766 

c 

100 

CONTINUE 

0767 

REWIND  9 

0768 

WRITE! 9)  ! ! AX ! J, I ),J=1, NBAND), 1=1, JQMAX) 

0769 

WRITEI9)  !X! I) ,1=1, JQMAX) 

0770 

END  FILE  9 

0771 

r 

REWIND  9 

0772 

L 

NBO  =  NBAND+1 

0773 

00  107  1=1,  JOMAX 

0774 

AX!NBD,I)  =  X!I) 

0775 

107 

CONTINUE 

C  WRITE(6,4000) 

C4000  FQRMATdHl,'  EQUATION',  35X,  'THE  BANDED  MATRIX  TERMS 
C  CALL  RITEd,  JQMAX,  NBO,  JOMAX,  NBD,  AX) 

C  WRITE(6,4003) 

C4003  FORMATdHl,  45X,  •***  THE  LOAD  VECTOR  Xd)  Iff  ) 

C  WRITE(6,4004)  (X(I),  1=1,  JQMAX) 

C4004  FORMAT(28(2X,  10D12.3  /  )) 

C 

CALL  TRMSTR(AX,  JQMAX,  NBD  ,  IBW,  IBW,  NBAND,  DT ,  RT 
C 

WRITE(6,4006)  ET,  RT ,  OT 

4006  FORMAT!///  '  ERROR  CONDITION  OF  SOLVER  ROUTINE  IS 
1  'RANK  IS  *,  F6.1,  5X,  'DETERMINANT  =  *,  G10.3) 
IF(ET.EQ,1.)  STOP  1 
C 

DO  108  1=1, JQMAX 
Xd)  =  AXd,I) 

108  CONTINUE 
C 

READ! 9)  (! AX! J,I) ,J=1, NBAND), 1=1, JQMAX) 

READ(9)  (R( I) ,1=1, JQMAX) 


OUTPUT  OF  THE  NODAL  DISPLACEMENTS,  U,  V,  W 


00010450 
00010460 
00010470 
00010480 
00010490 
00010500 
00010510 
00010520 
00010530 
00010540 
000105,50 
00010560 
00010570 
00010580 
00010590 
00010600 
00010610 
00010620 
00010630 
00010640 
00010650 
00010660 
00010670 
00010680 
00010690 
00010700 
00010710 
00010720 
00010730 
00010740 
,  AX! I,J) •  //)00010750 
00010760 
00010770 
00010780 
00010790 
00010800 
00010810 
,  ET)  00010820 

00010830 
00010840 
F4.1,  5X,  00010850 

00010860 
00010870 
00010880 
00010890 
00010900 
00010910 
,  00010920 

00010930 
00010940 
00010950 
00010960 
; « Jjc  *  >{c  ^  #  #  0  0  0 1  0  9  7  0 

00010980 
00010990 
00011000 
;#:«c:#:#ilcj)t:«t«»|tjjc3{cj!!#0001  1010 
00011020 
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0785 

0786 

0787 

0788 

0789 

0790 

0791 

0792 

0793 

0794 

0795 

0796 

0797 

0798 

0799 

0800 

0801 

0802 

0803 

0804 

0805 

0806 

0807 

0808 

0809 

0810 


0811 

0812 

0813 

0814 

0815 

0816 

0817 

0818 

0819 

0820 

0821 

0822 

0823 

0824 

0825 

0826 

0827 

0828 

0829 

0830 

0831 


C 

C 

c 

c 

c 

c 

c 


c 


c 


WRITe(6,650) 

J  =  1 

DO  12  IK  =  1,  LAW 
DO  11  JK  =  1,  LAT 

WRITE(6,651)  J,  XI3*J-2),  X(3*J~1),  XI3*J) 

J  =  J+1 

11  CONTINUE 
WRITE(6,653) 

12  CONTINUE 

WRITEt6t9950) 

9950  FORMATdHl,  5X,  ’EQUATION*,  5X,  ****  THE  ACCURACY  TEST, 
1  ****  ,  lOX,  ****  THE  AVERAGE  ABSOLUTE  ERROR  ///J 

ERR  =  0.00 
DO  9990  I=1,JQMAX 
TEST  =  0,00 
DO  9960  J=1,NBAND 
IC  -  I+J-IBW-1 
IF( I.LE.IBWl)  IC  -  J 
IF{ IC.GT.JQMAX)  GO  TO  9970 
TEST  =  TEST+AX(  J,n*X{  IC) 

9960  CONTINUE 

9970  TEST  =  TEST-R(I) 

ERR  =  ERR+DABS(TEST) 

AVE  =  ERR/I 

WRITE(6,9980)  I,  TEST,  AVE 
9980  F0RMAT(5X,  I4,10X,  G15.8,  32X,  G15.8) 

9990  CONTINUE 


00011030 
00011040 
00011050 
00011060 
00011070 
00011080 
00011090 
00011100 
OOOllIlO 
00011120 
00011130 
TEST-Rd)  00011140 
00011150 
00011160 
00011170 
00011180 
00011190 
00011200 
00011210 
00011220 
00011230 
00011240 
00011250 
00011260 
00011270 
00011280 
00011290 
00011300 


00011310 


CALCULATION  OF  THE  STRAIN  (S)  AND  STRESS  (T) 


00011330 

00011340 


00011350 

1360 


SXM  =  SXMAX  *  1.E06 
SXE  =  C3E  *  1.E06 
WRITE(6,670)  SXM,  SXE 
WRITE(6,671) 

HR  =  l./(2.*H) 

XK  =  FLOATIK) 


00011370 

00011380 

00011390 

00011400 

00011410 

00011420 

00011430 


DO  399  1=1,  LAW 
DO  398  J=l,  LAT 


00011440 

00011450 

00011460 


11=1-1 

12=1-2 

NODE  =  LAT’J'Il+J 
JJl  =  3*tLAT4Il+J)-2 
JJ2  =  3>«'(LAT*I2+J)-2 
JJ3  =  3*{LAT*I2+J)-5 
JJ4  =  3«(LAT*I+J)-2 
JJ5  =  3«aAT*I  +  J)  +  l 
JJ6  =  3*(LAT*I1+J)+1 
JJ7  =  3*(LAT*I2+J)+1 
JJ8  =  3#{ LAT*Il+J)-5 
JJ9  =  3*(LAT*I+J)-5 
JJIO  =3*(LAT>«'Il+J)-8 


00011470 

00011480 

00011490 

00011500 

00011510 

00011520 

00011530 

00011540 

00011550 

00011560 

00011570 

00011580 

00011590 

00011600 
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0832 

JJll  =3«(LAT*( I+l)+J)-2 

00011610 

0833 

JJ12  =3=0=(LAT*Il  +  J)+4 

00011620 

0834 

JJ13  =3*(LAT>«'(  I-3)+J)-2 

00011630 

C 

00011640 

0835 

Z  =  (FLOAT(  J)-(FL0AT(LAT)+1.  )/2. 

00011650 

0836 

SX  =  C2^Z  +  C3 

00011660 

C 

00011670 

0837 

IF( I.EQ.l)  GO  TO  385 

00011680 

0838 

IF( I. EQ. LAW)  GO  TO  386 

00011690 

0839 

IF(I.GT.FSW1.AND.I.LT.FSW2)  GO  TO  382 

00011700 

0840 

IF( I.EQ.FSWl.OR. I.EQ.FSW2)  GO  TO  383 

00011710 

C 

00011720 

0841 

HI  =  H 

00011730 

0842 

H2  =  HI 

00011740 

0843 

GO  TO  384 

00011750 

C 

00011760 

0844 

382 

HI  ==  XK«H 

00011770 

0845 

H2  =  HI 

00011780 

0846 

GO  TO  384 

00011790 

C 

00011800 

0847 

383 

HI  =  H 

00011810 

0848 

H2  =  XK^S'H 

00011820 

0849 

IF{ I.EO.FSWl)  GO  TO  384 

00011830 

0850 

HI  =  XK’i'H 

00011840 

0851 

H2  =  H 

00011850 

C 

00011860 

0852 

384 

H12  =  H1/H2 

00011870 

0853 

H21  =  H2/H1 

00011880 

0854 

HRD  =  (H2“H1)/{H1«H2) 

00011890 

0855 

HRS  =  l./(Hl+H2) 

00011900 

C 

00011910 

0856 

SY  =  HRS«IH12«X(  JJ4+1)-H21>5'X(  JJ2  +  1)  )+HRD’S‘X(  JJl+1)  +  DV’i^Z 

+  BV 

00011920 

0857 

SXY  =  HRS»>'(H12=«'X{  JJ4)-H21’»‘X{  JJ2)  )  +  HRD=!‘X(JJ1)  +  2.’5^C4=#=Z 

+  BU 

00011930 

0858 

SYZI  =  HRS=«'(H12=S'X(  JJ4+2)-H21«X(  JJ2+2)  )+HRD’J'X{  JJl  +  2) 

00011940 

0859 

GO  TO  387 

0001 1950 

C 

00011960 

0860 

385 

SY  =  HR*(4.>«'X(  JJ4+1}-3.*X(  JJ1+1)-X(  JJll  +  1)  )  +  DV«Z  +  BV 

00011970 

0861 

SXY  =  HR*(4.*X(  JJ4)-3.*X(  JJ1)-X(  JJll)  )  +  2.=!'C4’J'Z  +  BU 

00011980 

0862 

SYZI  =  HR*(4.=!=X(  JJ4+2)-3.*X(  JJl  +  2)-X(  JJll  +  2)  ) 

00011990 

0863 

GO  TO  387 

00012000 

C 

00012010 

0864 

386 

SY  =  HR^(3.#X( JJ1+1)+X( JJ13+1)-4.*X{ JJ2+1) )  +  DV*Z  +  BV 

00012020 

0865 

SXY  =  HR=«‘(3.=«'X(  JJ1)+X(  JJ13)-4.=<‘X{  JJ2)  )  +  4-  BU 

00012030 

0866 

SYZI  =  HR=i‘t3.=S'X(  JJl+2)+X(  JJ13  +  2)~4.*X{  JJ2  +  2)  ) 

00012040 

C 

00012050 

0867 

387 

DO  392  M=l,  NLAY 

00012060 

0868 

IF(M.EQ.l.AND.J.GT.INFIl))  GO  TO  392 

00012070 

0869 

IF(M.EQ.l)  GO  TO  388 

00012080 

0870 

IF(J.LE.INF(M-1).0R.J.GT.INF(M) )  GO  TO  392 

00012090 

0871 

388 

IF(J.EQ.l)  GO  TO  389 

00012100 

0872 

IF( J.EQ. INFJM) .OR. J.EQ.LAT)  GO  TO  390 

00012110 

C 

00012120 

0873 

SZ  =  HR4c(X(  JJ6+2)~X(  JJ8+2)  ) 

00012130 

0874 

SYZJ  =  HR»S'(  X(  JJ6+1  )-X(  JJ8+1  )  ) 

00012140 

0875 

SXZ  =  HR*( X( JJ6)-X{ JJ8)  ) 

00012150 

0876 

GO  TO  391 

00012160 

C 

00012170 

0877 

389 

SZ  =  HR=*'(4.>S‘X(  JJ6+2)-3.=i'X{JJl+2)-X(  JJ12  +  2)  ) 

00012180 
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0878 

0879 

0880 

0881 

0882 

0883 

0884 


0885 

0886 

0887 

0888 

0889 

0890 

0891 

0892 


0893 

0894 

0895 

0896 

0897 

0898 

0899 

0900 

0901 

0902 

0903 

0904 

0905 

0906 

0907 

0908 

0909 

0910 


0911 

0912 

0913 


SYZJ  =  HR*(4**X( JJ6+1)-3.«X{ JJ1+1)-X( JJ12+1) ) 

SXZ  =  HR*(4.*X(  JJ6)-3.*X{  JJ1)-XUJ12I  ) 

GO  TO  391 
C 

390  SZ  =  HR#(3,*X( JJl+2)+X(JJ10+2)-4.*X( JJ8+2) ) 

SYZJ  =  HR*(3.*X( JJ1+1)+X( JJ10+1)'4.*X( JJ8+1) } 

SXZ  =  HR*(3,*X{ JJ1)+X( JJ10)-4.*X(  JJ8> ) 

C 

391  SYZ  =  SYZI  +  SYZJ 
C 

C  CALCULATION  OF  THE  STRESS  (T) 

C 

TX  =  C11(M)>!'SX  +  C12{M)*SY  +  C13(M)*SZ  +  C16(M)>!'SXY 

TY  =  C12{M)*SX  +  C22(M)»SY  +  C23(M)«SZ  +  C26(M)’!=SXY 

TZ  =  C13(M)*SX  +  C23(M):«=SY  +  C33(M)’)'SZ  +  C36(M)^SXY 

C 

TYZ  =  C44(M)*SY2  +  C45{M)^SXZ 
TXZ  =  C45(M)*SYZ  +  C55(M)*SXZ 

^  TXY  =  C16{M)*SX  +  C26(M)=»'SY  +  C36(M)=i‘SZ  +  C66(M)=>'SXY 

WRITE(6,672)  NODE,  TX,  TY,  TZ,  TYZ,  TXZ,  TXY,  SY,  SZ,  SYZ, 
WRITE(6,397)  SX 
C 

C  STRESS  AND  STRAINS  JUST  ABOVE  AN  INTERFACE 
C 

IF( J.NE.  INF{M) .OR, J,EQ,LAT)  GO  TO  392 
P  =  M+1 

SZ  =  HR*(  4,’«'X(  JJ6+2  )-3.*X(  JJl+2)-X(  JJ12  +  2)  ) 

SYZJ  =  HR=^(4.=^X(  JJ6+1  )-3.*X  (  JJl  +  1  )-X(  JJ12  +  1  )  ) 

SXZ  =  HR«(4.«X{ JJ6}~3.#X( JJ1)-X( JJ12) ) 

SYZ  =  SYZI  +  SYZJ 
C 

TX  =  Cll(P)=frSX  +  C12(P)4*SY  +  C13(P)*SZ  +  C16(P)*SXY 

TY  -  C12(P)*SX  +  C22(P)#SY  +  C23{P)*SZ  +  C26(P)=<‘SXY 

tZ  =  C13(P)*SX  +  C23{P)’i‘SY  -f  C33{P)*SZ  +  C36{P)*SXY 

C 

TYZ  =  C44(P)«SYZ  +  C45(P)>«'SXZ 
TXZ  =  C45(P}<=SYZ  +  C55(P)*SXZ 

TXY  =  C16(P)^SX  +  C26(P)>;^SY  +  C36(P)*SZ  +  C66{P)=^SXY 
C 

^  WRITE(6,672)  NODE,  TX,  TY,  TZ,  TYZ,  TXZ,  TXY,  SY,  SZ,  SYZ, 

392  CONTINUE 
398  CONTINUE 

WRITE(6,652) 


00012190 

00012200 

00012210 

00012220 

00012230 

00012240 

00012250 

00012260 

00012270 

00012280 

00012290 

00012300 

00012310 

00012320 

00012330 

00012340 

00012350 

00012360 

00012370 

00012380 

SXZ,SXY00012390 

00012400 

00012410 

00012420 

00012430 

00012440 

00012450 

00012460 

00012470 

00012480 

00012490 

00012500 

00012510 

00012520 

00012530 

00012540 

00012550 

00012560 

00012570 

00012580 

SXZ,SXY00012590 

00012600 

00012610 

00012620 

00012630 


C 

C 

C 

C 

C 

C 

C 


C 


399  CONTINUE 
9000  CONTINUE 


00012640 

00012650 


00012660 

###  ^ ^  ^  5}c  5}c  «««««««««**  :ic  ##  j}, ^  0  0 0  1 2  6 7  0 


FORMATS 


00012680 

00012690 


* ♦  *  *  *  *  5}t  * J{c  >>c *  )J5  3{t  )(c  :fc « 3):  ;(i:  Jjt  Jjt  *  3{:  S}(  *  3}:  :(c  !>:  ^  Jit  jf:  *  jjs  J{:  3>;  sc  ;{£  :{£  s{:  *  *  Jje : 


00012700 
*  sC  :C  =>:  *  SC  sC  « sft  sC  sjt  sC 0 0 0  1  2  7  1  0 


397  FORMAT! 14X,1P1E11. 3/) 

600  FORMAT!  IHl,  44X,  44H*«*  UNIFORM  BENDING  OF  A  LAMINATED  PLATE  **>«=) 

601  FORMAT!5I10) 


00012720 

00012730 

00012740 

00012750 


00012760 
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0914 


0915 

0916 

0917 

0918 

0919 

0920 

0921 

0922 

0923 

0924 

0925 

0926 

0927 

0928 

0929 


0930 

0931 

0932 

0933 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


18X,  37HCHANGE  IN  MESH  WIDTH  IFSWl)  AT  I 
18X»  37HCHANGE  IN  MESH  WIDTH  (FSW2)  AT  I 


=  ,  14  // 
=  ,  14  // 


602  FORMAT(////  5X,  ISH^’^’S^  INPUT  DATA  ***  ///  00012770 

1  18X,  ‘NUMBER  OF  LAYERS  IN  CROSS  SECTIONr  NLAY  14  //  00012780 

2  18X,  ‘NUMBER  OF  NODES  ON  VERTICAL  AXIS»  LAT  14  //  00012790 

3  18X,  ‘NUMBER  OF  NODES  ON  HORIZONTAL  AXISt  LAW  -*♦  14  ///  00012800 

4  18X,  37HCHANGE  IN  MESH  WIDTH  IFSWl)  AT  I  =  ,  14  //  00012810 

5  18X»  37HCHANGE  IN  MESH  WIDTH  (FSW2)  AT  I  =  i  14  //  00012820 

6  18X,  37HMESH  WIDTH  MAGNIFICATION  FACTOR,  K  =  ,  14  /)  00012830 

00012840 

603  F0RMAT(8G12.5)  00012850 

00012860 

604  FORMATdHl,  55X,  21H**«  MATERIAL  DATA  ////  2X,  5HLAYER,  7X,  00012870 

1  3HE11,  9X,  3HE22,  9X,  3HE33,  9X,  3HE12,  9X,  3HE13,  9X,  00012880 

2  3HE23,  8X,  4HNU12,  4X,  4HNU13t  4X,  4HNU23  //  )  00012890 

00012900 

605  FGRMATOX,  12,  6X,  2PE10.3,  2(2X,  1PE10.3),  3{2X,  0PE10,3),  00012910 

1  3(3X,  F5.2)  /  )  00012920 

00012930 

606  FORMATUOG10.3)  00012940 

00012950 

607  FORMAT!///  18X,  26HFINE  SIMULATION  WIDTH,  H  ==  ,F8.5)  00012960 

00012970 

608  FORMAT!//  18X,  9HLAYER  NO.,  2X,  13,  5X,  17HINTERFACE  AT  J  ,13)  00012980 

00012990 

611  FORMAT!//  45X,  4lH>i'=S'*  COEFFICIENTS  OF  THERMAL  EXPANSION  ///  00013000 

1  IX,  5HLAYER,  8X,  5HTHETA,  12X,  3HALI,  i2X,  3HAL2,  12X,  00013010 

2  3HAL3,  12X,  3HAL6,  12X,  4HAL1P,  IIX,  4HAL2P,  IIX,  4HAL3P  00013020 

3  ///  )  00013030 

00013040 

613  FORMAT!//  53X,  26H*#*  STIFFNESS  MATRICES  ///  IX,  00013050 

1  IIHLAYER/THETA,  21X,  12HX-Y-Z  MATRIX,  44X,  00013060 

2  18HX-Y-Z  PRIME  MATRIX  ///  )  00013070 

00013080 

614  F0RMATI2X,  12,  9X,  F5.1,  5X,  7!5X,  E10.3))  00013090 

00013100 

620  F0RMAT(2X,  12,  5X,  1P12E10.3  //  19X,  5E10.3,  lOX,  5E10.3  //  29X,  00013110 

1  4E10.3,  20X,  4E10.3  //  1X,0PF5.1,  33X,  1P3E10,3,  30X,  00013120 

2  3E10.3  //  49X,  2E10.3,  40X,  2E10.3  //  59X,  ElO.3,  50X,  00013130 

3  E10.3  ///  )  00013140 

00013150 

650  FORMATIlHl  //  lOX,  GRID  POINT  DISPLACEMENT  FUNCTIONS  ///00013160 

1  16X,  5H  NODE,  5X,  14HU-DI SPLACEMENT,  6X,  14HV-D I SPLACEMENT , 00013170 

2  6X,  14HW-0ISPLACEMENT  ///  )  00013180 

00013190 

651  FORMATdOX,  110,  3E20.6  //  )  00013200 

652  FORMAT!//  12H  //  )  00013210 

653  FORMAT!//  lOX,  12H  //  )  00013220 

00013230 

670  FORMAT! IHl,  lOX,  77H**#  OUTPUT  STRESSES  AND  STRAINS  FOR  A  MAXIMUM  00013240 
ILONGITUDINAL  BENDING  STRAIN  OF  ,  F6.0,  22H  M ICRD-INCHES/INCH  AND  /00013250 

2  48X,  40H  AN  APPLIED  AXIAL  EXTENSIONAL  STRAIN  OF  ,  F6.0,  00013260 

3  19H  MICRO-INCHES/INCH.  //  lOX,  'NOTE:  INTERFACE  NODES  ARE  RE  PE ATOOOl 3270 

4ED  WITH  VALUES  GIVEN  BELOW  AND  ABOVE  THE  INTERFACE  RESPECTIVELY.'  00013280 
5  ////  )  00013290 

00013300 

671  F0RMATdX,5HN0DE  ,  5X,  5HSIG-X,  6X,  5HSIG-Y,  6X,  5HSIG-Z,  6X,  00013310 

1  6HTAU-YZ,  5X,  6HTAU-XZ,  5X,  6HTAU-XY,  5X,  SHEPS-Y,  6X,  00013320 

2  5HEPS-Z,  6X,  6HEPS-YZ,  5X,  6HEPS-XZ,  5X,  6HEPS-XY  /  17X,  00013330 

3  5HEPS-X  ///  )  00013340 

00013350 

-672  FORMATdX,  13,  6X,  1P11E11.3  /)  00013370 


6  18X,  37HMESH  WIDTH  MAGNIFICATION  FACTOR,  K  =  ,  14  /) 

603  F0RMAT!8G12.5) 

604  FORMATdHl,  55X,  21H**«  MATERIAL  DATA  ////  2X,  5HLAYER,  7X, 

1  3HE11,  9X,  3HE22,  9X,  3HE33,  9X,  3HE12,  9X,  3HE13,  9X, 

2  3HE23,  8X,  4HNU12,  4X,  4HNU13,  4X,  4HNU23  //  ) 

605  F0RMATI3X,  12,  6X,  2PE10.3,  2!2X,  1PE10.3),  3{2X,  0PE10.3), 

1  3(3X,  F5.2)  /  ) 

606  FORMATdOGlO.3) 

607  FORMAT!///  18X,  26HFINE  SIMULATION  WIDTH,  H  ==  ,F8.5) 

608  FORMAT!//  18X,  9HLAYER  NO.,  2X,  13,  5X,  17HINTERFACE  AT  J  ,13) 

611  FORMAT!//  45X,  41H>S'=S'*  COEFFICIENTS  OF  THERMAL  EXPANSION  /// 

1  IX,  5HLAYER,  8X,  5HTHETA,  12X,  3HAL1,  i2X,  3HAL2,  i2X, 

2  3HAL3,  12X,  3HAL6,  12X,  4HAL1P,  IIX,  4HAL2P,  IIX,  4HAL3P 

3  ///  ) 

613  FORMAT!//  53X,  26H*#*  STIFFNESS  MATRICES  ///  IX, 

1  llHLAYER/THETA,  21X,  12HX-Y-Z  MATRIX,  44X, 

2  18HX-Y-Z  PRIME  MATRIX  ///  ) 

614  F0RMATI2X,  12,  9X,  F5.1,  5X,  7!5X,  E10.3)) 

620  F0RMAT(2X,  12,  5X,  1P12E10.3  //  19X,  5E10.3,  lOX,  5E10.3  //  29X, 

1  4E10.3,  20X,  4E10.3  //  1X,0PF5.1,  33X,  1P3E10.3,  30X, 

2  3E10.3  //  49X,  2E10.3,  40X,  2E10.3  //  59X,  ElO.3,  50X, 

3  E10.3  ///  ) 


651  FORMATdOX,  110,  3E20.6  //  ) 

652  FORMAT!//  12H  //  ) 

653  FORMAT!//  lOX,  12H  //  ) 


STOP 

END 


00013380 

00013390 
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0001 

SUBROUTINE  MATCON 

00013400 

C 

00013410 

C 

j}c #  *  ijt #  «!c *  4: »!£ * «*  *  « « sic* # #  # «# *  « # #3)t *  # ##  # ;!c sS: *  # 5jc s}:  Jjc 5}: # # #  « :(c  # ##  # # * #  0  0 0  1  3  4  2 0 

C 

00013430 

C 

CALCULATION  OF  LAMINATE  LOAD  CONSTANTS  FOR  A  FULL  CROSS  SECTION 

00013440 

C 

00013450 

C 

**********slc5!c***sjt**sic*****!{c***sie***5}c**3}c********>)c*************************000  13460 

C 

00013470 

C 

THIS  SUBROUTINE  IS  GOOD  FOR  BENDING  OF  AN  ARBITRARILY  LAID  UP 

00013471 

C 

LAMINATE  WHICH  IS  SYMMETRIC  OR  NONSYMMETRIC  ABOUT  THE  MIDPLANE. 

00013472 

C 

00013480 

C 

THE  CONSTANTS  ARE  C2  =  INVERSE  BENDING  RADIUS 

00013490 

C 

C3E  =  APPLIED  UNIFORM  EXTENSIONAL  STRAIN 

00013500 

C 

C3  =  EXTENSIONAL  COUPLING  DUE  TO  BENDING  PLUS  C3E 

00013510 

C 

-  C4  =  IN-PLANE  SHEAR  COUPLING 

00013520 

C 

00013530 

C 

BU  OCCURS  IN  THE  -FGTN.  U{ Y,Z) 

00013540 

C 

BV  AND  DV  OCCUR  IN  THE  FCTN.  V(Y,Z) 

00013550 

C 

00013560 

C 

SXMAX  (EFFECTIVELY  THE  LOAD  INPUT)  IS  A  MAXIMUM  STRAIN 

00013570 

C 

00013580 

0002 

INTEGER  ORDER 

00013590 

c 

00013600 

0003 

COMMON  /MC/  ClK 6) » C12l 6) »C16( 6) tC22 ( 6 ) ,C26( 6) ,C66( 6) ,C13 (6 ) , 

00013610 

1  C23(6) ,C36(6) ,C44(6) tC45(6),C55{6)»C33{6) ,AL1(6) ,AL2(6) , 

00013620 

2  AL3(6) ,AL6(6) ,C2»C3,C3E,C4TBUtDU,BV,DV,HTSXMAX,NLAYTlNF(6) 

00013630 

c 

00013640 

0004 

DIMENSION  A(3,3),  B(3»3),  D(3,3)t  QM(3,3) 

00013650 

c 

00013660 

0005 

DOUBLE  PRECISION  A,  B,  D 

00013670 

c 

00013680 

0006 

ORDER  -  3 

00013690 

c 

00013700 

0007 

LAY  =  INF{1)-1 

00013710 

0008 

HL  =  H’i'FLOATtLAY) 

00013720 

0009 

HL2  =  HL«42/2. 

00013730 

0010 

HL3  =  HL’!'*3/3. 

00013740 

0011 

RN  =  FLOAT(NLAY) 

00013750 

0012 

RN2  =  RN*«2 

00013760 

c 

00013770 

0013 

DO  20  1=1,3 

00013780 

0014 

DO  20  J=l,3 

00013790 

0015 

A(I,J)  =  O.DO 

00013800 

0016 

B(  I,  J)  =  O.DO 

00013810 

0017 

D(I,J)  =  O.DO 

00013820 

00l8 

20  CONTINUE 

00013830 

c 

00013840 

0019 

DO  30  1=1,3 

00013850 

0020 

DO  30  J=l,3 

00013860 

0021 

DO  30  M=1,NLAY 

00013870 

0022 

QM(1,1)  =  C11(M)-C13(M)4C13(M)/C33(M) 

00013880 

0023 

QM(1,2)  =  C12(M)-C13{M)*C23(M)/C33(M} 

00013890 

0024 

QM(1,3)  =  C16(M)-C13(M)=S=C36{M)/C33(M) 

00013900 

0025 

QM{2,1)  =  QM(1,2) 

00013910 

0026 

0M(2,2)  =  C22(M)-C23(M)*C23(M)/C33(M) 

00013920 

0027 

QM{2,3)  =  C26(M)-C23(M)*C36(M)/C33(M) 

00013930 

0028 

QM(3,1)  =  QM(1,3) 

00013940 

0029 

0M(3,2)  =  QM(2,3) 

00013950 
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0030 

qM(3,3)  C66(M)-C36(M)=;.C36(M)/C33(M) 

00013960 

C 

OOOI397O 

C 

NOTE  THAT  THE  SUBSCRIPT  3  IN  QM  REPLACES  A  6  IN  STANDARD  NOTATION. 

00013980 

C 

THE  SAME  IS  TRUE  BELOW  IN  AllfJlt  B(I,J),  DII^J),  ETC. 

00013990 

C 

00014000 

0031 

Ml  =  2=J=M~1 

00014010 

0032 

M2  =  3’i‘M>S'(M-l  )  +  l 

00014020 

C 

00014030 

0033 

AU  ,  J)  =  A(  It  J)  +  HL^QMC  I,  J) 

00014040 

0034 

B(I,J)  =  BdtJ)  +  HL2=!‘QM{  I,  J>-(M1-RN) 

00014050 

0035 

DdtJ)  =  D(ItJ)  +  HL3-QM(  I,J)«(M2-1.5*RN*Ml  +  .75=i'RN2) 

00014060 

0036 

30. CONTINUE 

00014070 

C 

00014080 

C 

INVERT  (A).  STORE  IN  (A). 

00014090 

C 

00014100 

0037 

CALL  MATIN4  (AtORDER) 

00014110 

C 

00014120 

C 

MULTIPLY  (A)  INVERSE  ^  (B).  STORE  IN  A. 

00014130 

C 

00014140 

0038 

CALL  MAMULT  ( A » B t ORDER, A ) 

00014150 

C 

00014160 

C 

MULTIPLY  (B)  ’J'  (A)INVERSE  =5=  (B).  STORE  IN  B. 

00014170 

C 

00014180 

0039 

CALL  MAMULT  { B , A , ORDER , B ) 

00014190 

C 

00014200 

0040 

DO  40  1=1,3 

00014210 

0041 

DO  40  J=l,3 

00014220 

0042 

A(  I  ,J}  =  -l.’^AI  ItJ) 

00014230 

0043 

D( I,J)  =  D(  I,J)  -  B(  ItJ) 

00014240 

0044 

40  CONTINUE 

00014250 

C 

00014260 

C 

INVERT  NEW  MATRIX  (D).  THE  RESULT  IS  D-PRIME.  STORE  IN  D. 

00014270 

C 

00014280 

0045 

CALL  MATIN4  (DtORDER) 

00014290 

c 

00014300 

c 

MULTIPLY  -(A) INVERSE  ^  B  ^  D-PRIME  WHICH  YIELDS  B-PRIME.  STORE  IN  B. 

00014310 

c 

00014320 

0046 

CALL  MAMULT  ( A t D  ,ORDER , B ) 

00014330 

c 

00014340 

c 

DETERMINE  THE  LOAD  CONSTANTS.  MINUS  C2  IMPLIES  A  SMILING  PLATE. 

00014350 

c 

00014360 

0047 

ZMAX  =  RN*:'HL/2. 

00014370 

0048 

C2  =  -D(  It  1)=!=SXMAX/{B(  Itl)  +D  (  1 ,  1 )  =?=  ZMAX  ) 

00014380 

0049 

RATIO  =  C2/D{1,1) 

00014390 

c 

00014400 

0050 

C3  =  B(  1, 1  )=:'RATIO  +  C3E 

00014410 

0051 

C4  =  lt3)=!=RATI0 

00014420 

0052 

BU  =  B(3tl)-'=RATI0 

00014430 

0053 

BV  =  B(2tl)*RATI0 

00014440 

0054 

DV  =  0(  lt2)=4'RATI0 

00014450 

c 

RATIO  =  -RATIO 

00014460 

0055 

WRITE(6,50) 

00014470 

0056 

50  FORMAT!/////  48Xt  35H«=«=^  THE  LAMINATE  LOAD  CONSTANTS  ///  ) 

00014480 

0057 

WRITE{6,60)  C2t  C3t  C4,  BU,  BV,  DV,  RATIO 

00014490 

0058 

60  FORMAT!'  C2  =  IPEIO.3,  4X,  *C3  =  E10.3,  4X,  'C4  =  ',E10.3, 

00014500 

1  4Xt  ‘BU  =  ',  ElO.3,  4X,  'BV  =  E10.3,  4X,  * DV  =  *,  E10.3 

,00014510 

2  4Xt  'MT  =  • t  E10.3  ) 

00014520 

0059 

RETURN 

00014530 

0060 

END 

00014540 
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^OPTIONS  IN  EFFECT*  NOTERM , NO  ID , EBCD IC , SOURCE , NOL I  ST , NODECK , LOAD , NOMAP , NOTEST 
^OPTIONS  IN  EFFECT*  NAME  =  MATCON  »  LINECNT  =  60 

*STATISTICS*  source  STATEMENTS  =  60, PROGRAM  SIZE  =  2060 

*STATISTICS*  NO  DIAGNOSTICS  GENERATED 
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OOOl 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 


SUBROUTINE  MAMULT ( B , C , N , A ) 

C 

C  MAMULT  POSTMULTI PL lES  MATRIX  (B)  BY  MATRIX  (C)  AND  STORES  THE 
C  RESULT  IN  MATRIX  (A)  WHERE  N  IS  THE  ORDER  OF  THE  MATRICES. 

C 

DOUBLE  PRECISION  A,B,C,SUM 
DIMENSION  A(N,N),  6(N,N),  C(N,N) 

DO  1  1=1, N 
DO  1  J=1  ,N 
SUM  =  0. 

DO  2  K=1,N 

SUM  =  SUM  +  B(  I,K}*C(K, J} 

2  CONTINUE 
A( I, J)  =  SUM 
I  CONTINUE 
RETURN 
END 


00014550 

00014551 

00014552 

00014553 

00014554 

00014560 

00014570 

00014580 

00014590 

00014600 

00014610 

00014620 

00014630 

00014640 

00014650 

00014660 

00014670 
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*OPTIONS  IN  EFFECT*  NOTERM, NOI D , EBCD IC , SOURCE , NOL I  ST, NODECK , LOAD , NDMAP , NOTEST 
^OPTIONS  IN  EFFECT*  NAME  =  MAMULT  ,  LINECNT  =  60 

*STATISTICS*  SOURCE  STATEMENTS  =  13, PROGRAM  SIZE  =  702 

*STATISTICS*  NO  DIAGNOSTICS  GENERATED 
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OOOl 

SUBROUTINE  MAT I N4 ( ARRAY , N ) 

00014680 

C 

00014681 

C 

MATIN4  INVERTS  THE  MATRIX  (ARRAY)  WHICH  IS  OF  ORDER  N. 

00014682 

C 

00014683 

0002 

DIMENSION  ARRAY(N,N) 

00014690 

0003 

DOUBLE  PRECISION  ARRAY 

00014700 

0004 

DO  604  1=1, N 

00014710 

0005 

STORE  =  ARRAY! 1,1) 

00014720 

0006 

ARRAY! 1,1)  =  1, 

00014730 

0007 

DO  601  J=1,N 

00014740 

0008 

601 

ARRAY(I,J)  =  ARRAY! I,J)/STORE 

00014750 

0009 

DO  604  K=1,N 

00014760 

0010 

IF(K-1)602,604,602 

00014770 

0011 

602 

STORE  =  ARRAYIK, I ) 

00014780 

0012 

ARRAYIK,  I  )  i:  0. 

00014790 

0013 

DO  603  J=1,N 

00014800 

0014 

603 

ARRAY(K,J)  =  ARRAY(K,J)  -  STORE*ARRAY ! I , J ) 

00014810 

0015 

604 

CONTINUE 

00014820 

0016 

RETURN 

00014830 

0017 

END 

00014840 
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0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 


SUBROUTINE  TRMSTR ( A,N,ND,NLDtNRD,NED,D,R,E ) 

C 

C  TRMSTR  IS  THE  SUBROUTINE  TRIMSS  WITH  MATRIX  A  TRANSPOSED. 

C  THE  SIMULTANEOUS  SOLUTIONS  IS  GAUSSIAN  ELIMINATION, 

C  MODIFIED  TO  TAKE  ADVANTAGE  OF  THE  REDUCED  MATRIX.  THE 

C  ROUTINE  ALSO  USES  PARTIAL  PIVOTING  TO  REDUCE  ROUNDOFF  ERROR. 

C  INPUT 

C  1  A 

C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 

c  2  N 

C  3  ND 

C 

C  4  NLD 

C 

c 

c  5  NRD 

C 
C 

C  6  NED 

C  OUTPUT 

C  1  A 

C 

C  2D 

C  3  R 

C  4  £ 

C 
C 
C 
C 

C  SUBROUTINE  TRMSTR! A, N,ND, NLD, NRD, NED, D,R,E) 

DIMENSION  A(ND,1) 

DOUBLE  PRECISION  A,D,Y,W,S 
XI  =  1, 

LI  =  1 
E  =  0. 

R  -  0. 

D=l. 

ND1=NED+1 
M  =  NLD 
NM1=N-1 
DO  1  1=1, NMl 
IF( I.GT. (N-NLD) )M=M-1 
NN=I+M-1 
DO  2  II=I,NN 


FIRST  LOCATION  OF  COEFFICIENT  MATRIX, I. E.  A(l,l). 
THE  BAND  ELEMENTS  IN  EACH  ROW  MUST  BE  LEFT 
JUSTIFIED  AND  EXTEND  TO  THE  RIGHT  M  PLACES 
(M=MIN(N,NLD+NRD+1).  IF  IN  ANY  PARTICULAR  ROW 
THERE  ARE  ONLY  K  BAND  ELEMENTS  AND  K  IS  LESS 
THAN  M,  THEN  THE  M-K  RIGHT  MOST  ELEMENTS  OF  THAT 
ROW  WILL  BE  SET  TO  ZERO.  THE  ROW  WHOSE  LEFT 
MOST  COLUMN  IN  THE  FULL  BLOWN  MATRIX  CONTAINS 
A  NON^ZERO  ELEMENT  MUST  BE  THE  FIRST  ROW  OF  THE 
REDUCED  MATRIX  AND  ETC,  THE  COLUMN  TO  THE 
IMMEDIATE  RIGHT  OF  THE'  REDUCED  MATRIX  (FORMED  AS 
ABOVE)  MUST  CONTAIN  THE  RIGHT  HAND  SIDE  OF  THE 
EQUATION  SET  IN  QUESTION.  IT  SHOULD  NOW  BE 
OBVIOUS  THAT  AN  N  X  N+1  FULL  BLOWN  SYSTEM  WOULD 
BE  REDUCED  BY  THE  ABOVE  METHOD  TO  AN  N  X  M+1 
SYSTEM. 

NUMBER  OF  SIMULTANEOUS  EQUATIONS  TO  BE  SOLVED. 
VARIABLE  DIMENSION  INTEGER.  MUST  BE  EQUAL  TO 
ROW  DIMENSION  OF  A  IN  CALLING  PROGRAM. 

MAXIMUM  NUMBER  OF  BAND  ELEMENTS  T-O  THE  LEFT 
OF  PRINCIPAL  DIAGONAL  IN  ANY  ROW  OF  SYSTEM  TO 
BE  DETERMINED. 

MAXIMUM  NUMBER  OF  BAND  ELEMENTS  TO  THE  RIGHT 
OF  PRINCIPAL  DIAGONAL  IN  ANY  ROW  OF  SYSTEM  TO 
BE  DETERMINED. 

NED=MIN(N,NLD+NRD+1 ) 

THE  FIRST  COLUMN  OF  A  CONTAINS  THE  SOLUTION 
VECTOR. 

CONTAINS  DETERMINANT  OF  A. 

CONTAINS  RANK  OF  A. 

E=0,,  SOLUTION  O.K.  E=l.,  A  SINGULAR. 

E=2.,  SOLUTION  ATTEMPTED,  BUT  A  ILL  CONDITIONED 
OR  SINGULAR.  IN  THIS  CASE  SOLUTIONS  SHOULD  BE 
CHECKED  TO  ASSURE  VALIDITY. 


00014850 

00014860 

00014870 

00014880 

00014890 

00014900 

00014910 

00014920 

00014930 

00014940 

00014950 

00014960 

00014970 

00014980 

00014990 

00015000 

00015010 

00015020 

00015030 

00015040 

00015050 

00015060 

00015070 

00015080 

00015090 

00015100 

00015110 

00015120 

00015130 

00015140 

00015150 

00015160 

00015170 

00015180 

00015190 

00015200 

00015210 

00015220 

00015230 

00015240 

00015250 

00015260 

00015270 

00015280 

00015290 

00015300 

00015310 

00015320 

00015330 

00015340 

00015350 

00015360 

00015370 

00015380 

00015390 

00015400 

00015410 

00015420 
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0016 

IF(DABS(A{ 1,1) ) •GE.DABS(A (1, I Ul) ) )  GO  TO  2 

00015430 

0017 

D=-D 

00015440 

0018 

DO  3  J=1,ND1 

00015450 

0019 

Y=A( J, I) 

0001546*0 

0020 

A( J, I  )=A( J,II  +  1) 

00015470 

0021 

3 

A( J,  II  +  1)=Y 

00015480 

0022 

2 

CONTINUE 

00015490 

C 

D=D*A( 1, I ) 

00015500 

0023 

IF(A(1,I)  .EQ.  0.)  GO  TO  10 

00015510 

0024 

GO  TO  (5, 13), LI 

00015520 

0025 

13 

IF (DABS (DABS ( (X1-A(1,I) ) /XI ) -1. ) • LT • 1 , E-07 ) 

E  =  2. 

00015530 

0026 

XI  =  A(1,I) 

00015540 

0027 

5 

R  =  R  +  1. 

00015550 

0028 

LI  =  2 

00015560 

0029 

DO  4  J=2,ND1 

00015570 

0030 

4 

A(J,I)=A(J,I  )/  A(1,I) 

00015580 

0031 

K=I  +  1 

00015590 

0032 

NN=I+M 

00015600 

0033 

DO  1  II=K,NN 

00015610 

0034 

W=A(1, II) 

00015620 

0035 

DO  6  J=1,NED 

00015630 

0036 

6 

A(J,II )=A(J+1,II)-A(J+1,I)*W 

00015640 

0037 

A(ND1,II)=A(NE0,II) 

00015650 

0038 

1 

A{NED, 1 1 )=0. 

00015660 

0039 

IF( A( 1,N) .EQ.O. )G0  TO  10 

00015670 

0040 

IF(DABS(DABS( 1 Xl-A( 1,N) ) /XI ) -1 . ) . LT, 1 . E>07 ) 

E=2. 

00015680 

0041 

9 

R  =  R  +  1, 

00015690 

0042 

A( 1,N}=A(ND1,N)/A( 1,N) 

00015700 

0043 

K  =  NM1 

00015710 

0044 

NN=2 

00015720 

0045 

8 

IF(NN,GT.NED)NN=NED 

00015730 

0046 

J  =  K+l 

00015740 

0047 

S  =  0. 

00015750 

0048 

DO  7  I=2,NN 

00015760 

0049 

S=S+A(1,J)*A( I,K) 

00015770 

0050 

7 

J  =  J+1 

00015780 

0051 

A(1,K)=AIND1,K)-S 

00015790 

0052 

NN=NN+1 

00015800 

0053 

K=K-1 

00015810 

0054 

IF(K,NE.O)GD  TO  8 

00015820 

0055 

RETURN 

00015830 

0056 

10 

E=l, 

00015840 

0057 

RETURN 

00015850 

0058 

END 

00015860 
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0001 

SUBROUTINE  R I TE ( I OUM, NR, NC , MR, MC , A ) 

00015870 

0002 

DOUBLE  PRECISION  A 

00015880 

0003 

DIMENSION  A(MR,MC) 

00015890 

0004 

IPRINT=  12 

00015900 

0005 

IF(  lOUM.NE.l)  IPRINT=  30 

00015910 

0006 

IPR=  IPRINT-1 

00015920 

0007 

DO  35  K=1,NC,IPRINT 

00015930 

0008 

MAX=  K+IPR 

00015940 

0009 

IF(MAX.GT.NC)  MAX=NC 

00015950 

0010 

IF(K.NE.l)  WRITE(6,103) 

00015960 

0011 

45 

WRITE(6,102)  {I,I=K,MAX) 

00015970 

0012 

DO  40  J=1,NR 

00015980 

0013 

40 

WRITE{6,105)  J,(A( Jtl) ,I=K,MAX) 

00015990 

0014 

35 

CONTINUE 

00016000 

0015 

RETURN 

00016010 

0016 

101 

F0RMAT(6X,30I4) 

00016020 

0017 

102 

F0RMAT(6X, 12110) 

00016030 

0018 

103 

FORMAT! »1» ) 

00016040 

0019 

104 

FORMAT!*  ',15,3014) 

00016050 

0020 

105 

FORMAT!'  *,15,12010.3) 

00016060 

0021 

END 

00016070 
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